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Abstract 

Ferroelectric  materials,  such  as  PZT,  PLZT  and  BaTi03,  are  being  considered,  or  are  already 
being  employed,  for  a  large  number  of  applications  including  nanopositioning,  high  speed  valves  for 
fuel  injectors,  ultrasonic  transducers,  high  speed  camera  shutters  and  auto  focusing  mechanisms, 
energy  harvesting,  and  pico  air  vehicle  design.  Their  advantages  include  nanometer  positioning 
resolution,  broadband  frequency  responses,  moderate  power  requirements,  the  capability  for  minia¬ 
turization,  and  complementary  actuator  and  sensor  capabilities.  However,  they  also  exhibit  creep, 
rate-dependent  hysteresis,  and  constitutive  nonlinearities  at  essentially  all  drive  levels  due  to  their 
noncentrosymmetric  nature.  In  this  paper,  we  model  the  hysteretic  dependence  of  strains  and  polar¬ 
ization  on  input  fields  and  stresses  using  the  homogenized  energy  model  (HEM)  framework.  At  the 
domain  level,  the  minimization  of  Gibbs  energy  densities  yields  linear  constitutive  relations.  Non- 
linearities  and  hysteresis  due  to  dipole  switching  is  modeled  at  the  grain  level  by  using  Boltzmann 
theory  to  specify  the  evolution  of  dipole  fractions  which  serve  as  internal  variables.  In  the  final  step 
of  the  development,  stochastic  homogenization,  based  on  the  assumption  that  interaction  fields  and 
driving  forces  are  manifestations  of  underlying  densities,  is  used  to  construct  nonlinear  constitutive 
relations  for  the  bulk  material.  It  is  demonstrated  that  these  relations  are  amenable  to  subsequent 
development  of  distributed  system  models.  The  paper  includes  significant  discussion  regarding  the 
mechanisms  that  produce  hysteresis  in  ferroelectric  materials.  The  capability  of  the  framework  for 
characterizing  various  hysteretic  phenomena,  including  creep  and  various  rate-dependencies,  is  illus¬ 
trated  by  validation  with  PZT  and  PLZT  data. 
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Nomenclature 


d  Piezoelectric  constant  (m/V  =  C/N) 

h  Piezoelectric  constant  (V/m  =  N/C) 

E  Electric  field  (V / m) 

g  Gibbs  energy  for  dipoles  (CV) 

Ga  Gibbs  energy  density  of  a- variant  (CV/rn3) 

Gea  Electric  Gibbs  energy  density  of  a- variant  (CV/rn3) 

P  Polarization  (C/m2) 

P  Polarization  kernel  (C/m2) 

P°  Polarization  of  a- variant  (C/m2) 

Pr  Remanent  polarization  of  a- variant  (C/m2) 

P ^  Minimum  polarization  of  a- variant  (C/m2) 

Elastic  compliance  of  a-variant  at  constant  field  (m2/N) 

T  Temperature  (K) 

x+,  x_,  X90  Fraction  of  positively,  negatively  and  90°  oriented  dipoles  (Unitless) 
yop  Elastic  stiffness  of  a-variant  at  constant  polarization  (N/m2) 

e  Permittivity  (F/m  =  C/Vrn) 

e  Strain  (Unitless) 

e  Strain  kernel 

ea  Strain  of  a-variant 

ef  Remanence  strain  of  a-variant 

Minimum  strain  of  a-variant 
7  7  =  V/  kT 

if  Inverse  susceptibility  at  constant  strain  (m/F  =  Vrn/C) 

c 7  Stress  (N/m2) 

790,  Tiso  Relaxation  time  for  90°  and  180°  switching  (s) 

Xe  Electric  susceptibility  (Unitless) 

Xa  Ferroelectric  susceptibility  of  a-variant  at  constant  stress  (F/m  =  C/Vrn) 

f>a  Helmholtz  energy  density  of  a-variant  (CV/rn3) 

1  Introduction 

Piezoelectric  materials  exhibit  two  complementary  properties  due  to  their  noncentrosymmetric  struc¬ 
ture:  the  direct  effect  in  which  applied  stresses  generate  an  electric  charge,  and  the  converse  effect 
in  which  applied  fields  produce  deformations  in  the  material.  These  properties  respectively  imbue 
the  materials  with  sensor  and  actuator  capabilities  as  well  as  the  potential  for  self-sensing  actuation. 
Commonly  employed  piezoelectric  materials  include  barium  titanate  (BaTiCU),  lead  zirconate  ti- 
tanate  or  PZT  (Pb(Tii__rZx)03)  lanthanum-doped  lead  zirconate  titanate  (PLZT),  and  piezoelectric 
polymers  such  as  polyvinylidene  fluoride  (PVDF).  For  the  PZT  compounds,  x  £  [0, 1]  is  chosen  to  op¬ 
timize  electromechanical  coupling.  Additionally,  naturally  occurring  crystals  such  as  quartz,  sucrose 
(table  sugar),  and  Rochelle  salts,  and  biological  materials  such  as  bone  exhibit  piezoelectric  effects 
to  varying  degrees.  The  actuator  and  sensor  capabilities  of  quartz  are  exemplified  by  the  fact  that  in 
1916,  Paul  Langevin  developed  a  quartz-based  sonar  transducer  for  submarine  detection  [64,75].  For 
present  applications,  however,  we  will  focus  on  BaTiOs,  PZT,  PLZT,  PVDF  and  the  electrostrictive 
material  lead  manganese  niobate  (PMN). 
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These  materials  are  ferroelectric  and  hence  exhibit  a  domain  structure  and  spontaneous  polariza¬ 
tion  at  temperatures  below  the  Curie  point.  The  metastability  associated  with  multiple  stable  dipole 
configurations  produces  hysteresis  and  constitutive  nonlinearities  in  field-polarization  and  field-strain 
responses.  In  this  paper,  we  focus  on  mechanisms  that  produce  hysteresis  in  ferroelectric  materials 
and  the  development  of  homogenized  energy  models  that  can  be  used  for  material  characterization, 
device  design,  and  model-based  control  design. 

There  exist  a  wide  range  of  actuator  and  sensor  designs  with  specific  choices  dictated  by  the 
application.  Stacked  actuators,  such  as  that  depicted  in  Figure  1(a),  are  employed  in  numerous 
applications  including  stages  for  atomic  force  microscopy  and  nanopositioners  (see  Figure  2),  high 
speed  valves  for  fuel  injection,  vibration  control  devices,  depth  finders  and  hydrophones,  and  linear 
and  rotary  piezomotors.  Tube  actuators  comprise  the  active  mechanisms  in  micropumps  and  scanning 
and  atomic  force  microscopes.  Bender-type  transducers,  such  as  the  unimorph  and  bimorph  designs 
depicted  in  Figure  1(c)  and  (d),  are  employed  in  pneumatic  values,  high  speed  camera  shutters, 


Interdigitated  Electrodes 


Figure  1:  (a)  Stack,  (b)  tube,  (c)  unimorph,  and  (d)  bimorph  transducers,  (e)  Linear  and  (f)  rotary 
piezomotors,  (g)  Macro-fiber  composite  (MFC)  transducer. 
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Figure  2:  (a)  Nanopositioning  stage,  (b)  high  speed  valve  for  fuel  injection,  (c)  inkjet  printer  nozzle, 
(d)  ultrasonic  transducer,  and  (e)  MFC  for  shape  modification  and  flow  control. 

energy  harvesting  devices,  piezoelectric  transformers,  and  inkjet  printers  (Figure  2(c)).  Bimorphs 
are  also  employed  as  actuators  for  pico  air  vehicles  such  as  the  Harvard  RoboBee  [90]  and  both 
flying  and  ambulatory  microrobots  [28,68,91].  Bending-type  PZT  actuators  are  also  employed  in 
ultrasonic  transducers  for  dental  tools  and  biomedical  imaging  and  treatment  (Figure  2(d)).  An 
emerging  technology  is  MRI-guided  focused  ultrasound  surgery  in  which  high  energy  ultrasound 
waves  are  used  to  thermally  ablate  tumors  and  growths  such  as  uterine  fibroids.  We  note  that 
ultrasonic  devices  operate  at  around  40  kHz  whereas  piezoelectric  transformers  operate  in  the  100  kHz 
to  1  MHz  range  thus  demonstrating  the  high  frequency  capabilities  of  the  materials.  Solid  state 
piezomotors  are  employed  in  applications  such  as  camera  auto-focusing  mechanisms  and  medical 
equipment  subject  to  large  magnetic  fields;  e.g.,  small  piezomotors  for  microsurgery  and  large  motors 
to  position  positions  in  MRI  environments.  Macro-fiber  composites  (MFC),  depicted  in  Figure  1(g), 
provide  large  strain  and  force  capabilities  in  addition  to  durability  and  flexibility.  They  are  presently 
being  considered  for  applications  including  shape  modification  and  flow  control  for  micro-air  vehicles 
(MAV)  as  well  as  energy  harvesting  in  a  range  of  environments  [12,49,82].  Sensor  applications  include 
accelerometers,  knock  sensors  to  monitor  engine  combustion,  pressure  and  force  sensors,  ultrasonic 
distance  sensors,  and  vibration  sensors  to  monitor  automotive,  railroad,  and  aircraft  components. 
The  direct  piezoelectric  effect  is  also  utilized  in  high  voltage  spark  igniters.  Finally,  180°  switching 
in  ferroelectric  materials  forms  the  basis  for  ferroelectric  memory  technologies  (e.g.,  FeRAM)  [70]. 
Additional  discussion  of  applications  can  be  found  in  [53,67,75,85-87]. 

The  advantages  of  piezoelectric  actuators,  sensors,  and  motors  are  due  to  a  number  of  factors. 
They  can  be  designed  to  provide  nanometer  positioning  resolution  and  operate  at  frequencies  ranging 
from  DC  to  MHz.  They  have  modest  power  requirements  and  do  not  create  nor  are  they  influenced 
by  electromagnetic  interference.  Their  solid  state  nature  promotes  miniaturization  and  simplified 
designs  which  improves  product  performance  at  reduced  costs.  They  generate  little  heat,  are  non¬ 
flammable,  and  can  be  operated  in  a  vacuum.  Finally,  the  complementary  direct  and  converse 
effects  provide  the  materials  and  devices  with  multiple  design  properties  including  actuator,  sensor, 
self-monitoring,  nondestructive  evaluation,  and  energy  harvesting  capabilities. 
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The  price  paid  for  the  unique  transducer  capabilities  of  ferroelectric  materials  is  hysteresis  and 
constitutive  nonlinearities  due  to  the  metastable  behavior  associated  with  inherent  domain  properties. 
Whereas  these  effects  can  be  minimized  using  charge  or  current  control  [62] ,  or  restriction  to  low  drive 
regimes,  this  can  increase  costs  and  limit  the  unique  transduction  capabilities  of  the  materials.  This 
necessitates  the  development  of  models  and  model-based  control  designs  that  quantify  the  nonlinear 
and  hysteretic  material  behavior  to  achieve  the  full  potential  of  the  materials  and  devices. 

Static,  quasistatic,  and  dynamic  hysteresis  behavior  that  must  be  incorporated  in  models  is 
illustrated  in  Figure  3.  The  field-polarization  and  field-strain  data  from  [96]  illustrates  that  rate- 
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Figure  3:  Rate-dependent  (a)  polarization  and  (b)  strain  PZT  data  from  [96].  (c)  Field-polarization, 
(d)  field-strain  and  (e)  time-strain  PZT  data  from  [92].  (f)  Field-strain  MFC  data  from  [32], 
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dependent  effects  are  significant  at  frequencies  as  low  as  1  Hz.  PZT  data  from  [92]  illustrates  that 
for  fixed  field  inputs,  both  the  strain  and  polarization  exhibit  significant  creep  on  timescales  of 
1  to  20  seconds.  Finally,  the  MFC  data  from  [32]  illustrates  nested  minor  loop  behavior  typical  of 
moderate  drive  regimes.  Whereas  the  full  switching  behavior  shown  in  Figure  3(a)-(d)  will  typically 
not  be  encountered  in  applications,  general  models  must  account  for  the  full  range  of  rate-dependent 
and  creep  behavior  to  provide  comprehensive  device  characterization. 

For  regimes  in  which  stress  and  field  inputs  a  and  E  are  maintained  at  low  levels,  the  nonlinear 
and  hysteretic  behavior  of  the  strain  e  and  polarization  P  can  be  approximated  by  the  linear  relations 

P  =  da  +  xaE 

E  (D 

£  =  sha  +  cLE 

originally  developed  by  Voigt  [88].  Here  d,xa  and  sE  denote  the  piezoelectric  constant,  ferroelec¬ 
tric  susceptibility  at  constant  stress,  and  elastic  compliance  at  constant  field  E.  These  are  often 
termed  the  piezoelectric  relations  which  infers  a  linear  connotation  on  piezoelectric  materials.  This 
is  somewhat  of  a  misnomer  since  the  piezoelectric  materials  BaTi03,  PZT,  PLZT,  and  PVDF  are 
ferroelectric  and  hence  exhibit  hysteresis  and  constitutive  nonlinearities.  For  some  applications, 
however,  the  linear  relations  (1)  provide  sufficient  accuracy,  especially  when  combined  with  feedback 
algorithms. 

To  motivate  various  nonlinear  modeling  strategies,  we  consider  the  multiscale  depiction  of  a  PZT- 
based  MFC  tranducer  in  Figure  4.  The  largest  scale  is  comprised  of  the  application  whereas  initial 
characterization  experiments  are  often  conducted  with  individual  MFC  bonded  to  an  elastic  structure; 


Figure  4:  Multiscale  behavior  of  a  PZT-based  MFC  transducer  at  the  application,  device,  material 
and  unit  cell  levels. 
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e.g.,  a  beam,  plate  or  shell.  The  active  component  of  the  MFC  is  PZT  fibers  which  are  depicted 
at  the  material,  grain,  domain  and  unit  cell  levels.  The  goal  for  device  optimization  and  control 
is  to  develop  nonlinear,  macroscale  constitutive  relations  that  can  be  employed  in  mechanical  and 
electrostatic  field  relations  which  subsequently  can  be  used  to  construct  finite  element  simulation  and 
control  codes  for  the  structure,  device  or  application.  The  different  modeling  hierarchies  are  defined 
by  the  degree  to  which  molecular,  domain,  and  grain-level  material  behavior  is  used  to  construct  the 
constitutive  relations. 

The  general  form  of  constitutive  relations  can  be  motivated  by  the  ionic  behavior  of  the  unit  cell. 
For  small  field  or  stress  inputs,  ionic  displacements  are  reversible  and  approximately  linear  so  they 
have  the  form 

po  =  da(J  +  X«E  +  pa 

£°  =  sEc  +  daE  +  ejf 

where  a  designates  the  dipole  variant  —  e.g.,  ±180,90  for  tetragonal  materials  —  and  are 

remanent  values  of  the  polarization  and  strain.  At  larger  input  levels,  however,  irreversible  switching 
occurs  which  produces  hysteresis  and  constitutive  nonlinearities.  This  combination  of  reversible  and 
irreversible  behavior  motivates  the  general  constitutive  formulation 

P(E ,  a)  =  d(E,  a)a  +  xaE  +  Pirr{E,  a) 

„  (3) 

e(E,  a)  =  sE a  ±  d(E,  a)E  ±  £irr(E,  a) 

where  d(E,  a),  Pjrr(E,  a)  and  £irr(E,  a)  incorporate  the  nonlinear  and  irreversible  history-dependence 
due  to  dipole  switching.  We  note  that  many  authors  employ  the  notation  Pr,£r  and  terminology 
‘irreversible  remanent  behavior’  but  the  concepts  are  the  same. 

Modeling  hierarchies  can  generally  be  defined  by  the  manner  in  which  d,  Pirr  and  £irr  are  con¬ 
structed.  Micromechanical,  or  microscopically-motivated  models  are  based  on  an  energy  description 
of  the  material  at  the  domain  or  grain  level  in  combination  with  various  homogenization  techniques 
to  provide  expressions  for  the  nonlinear,  effective  components  d,  Pirr  and  £irr.  The  objective  is  to  let 
the  micro-level  physics  inform,  to  the  degree  possible,  resulting  macroscale  behavior.  The  difficulty 
with  this  approach  is  that  resulting  models  are  often  too  complex  for  applications  requiring  high 
speed  implementation;  recall  that  applications  can  occur  at  kHz  rates.  Phenomenological  models 
circumvent  the  difficulties  associated  with  quantifying  complex,  or  poorly  understood,  micro-scale 
physics  by  constructing  relations  for  d,  Pirr  and  £irr  based  on  macroscale  observations  or  experi¬ 
mental  measurements.  In  many  models,  the  derivation  of  these  effective  components  is  guided  by 
thermodynamic  constraints,  and  the  history-dependence  associated  with  dipole  switching  is  often 
incorporated  with  nonphysical  internal  variables.  In  general,  phenomenological  models  are  less  com¬ 
plex  than  nricronrechanical  models  and  hence  are  amenable  to  finite  element  implementation  for 
devices  and  applications. 

Micromechanical  Models 

As  illustrated  in  Figure  4,  microscopically-motivated  or  nricronrechanical  models  must  incorporate 
two  phenomena:  a  kinematic  description  of  the  grain-level  irreversible  or  remanent  polarization 
Pirr  and  strain  £jrr,and  criteria  that  quantify  dipole  switching.  As  detailed  in  the  survey  papers 
[33,  51]  by  Huber  and  Landis,  many  early  papers  [36,  60]  assumed  single  domain  behavior  within 
single  crystals,  in  which  case,  the  remanent  polarization  and  strain  were  simply  the  spontaneous 
polarization  and  strain.  To  account  for  multiple  domains,  which  is  always  the  case  for  ferroelectric 
materials,  one  must  either  track  the  nucleation  and  movement  of  domain  walls,  or  account  for  their 
behavior  using  internal  variables  xa,  the  most  natural  of  which  are  volume  or  dipole  fractions. 
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Phase  field  models,  typically  based  on  Ginzburg-Landau  energy  relations,  provide  one  method  for 
quantifying  domain  wall  behavior  [89].  As  detailed  in  [75,81,95],  the  Ginzburg-Landau  theory  differs 
from  the  standard  Landau  relations  through  the  inclusion  of  polarization  gradient  energy  terms 
which  incorporate  local  polarization  changes  associated  with  domain  walls.  This  approach  permits 
fundamental  understanding  of  domain  behavior  but  is  typically  too  computationally  intensive  for  use 
in  macroscopic  constitutive  relations  requiring  high  speed  implementation.  To  reduce  the  number  of 
parameters,  continuum  theories  based  on  internal  variables  that  quantify  the  polarization  and  strain 
state  have  been  employed  by  a  number  of  researchers  —  see  [35]  as  well  as  the  summary  in  [33,51]. 

A  variety  of  techniques  have  been  employed  to  quantify  dipole  switching.  Hwang  et  al.  [36] 
introduced  the  concept  of  switching  levels  and  a  number  of  researchers  have  employed  switching 
surfaces  [73].  The  homogenized  energy  model  (HEM)  framework,  initially  developed  in  [75,78,80] 
and  extended  here,  quantifies  switching  through  the  evaluation  of  likelihoods  that  balance  the  Gibbs 
and  relative  thermal  energy  at  the  grain  level. 

We  note  that  the  kinetics  associated  with  dipole  switching  determine  the  time  scale  for  polar¬ 
ization  and  strain  changes  relative  to  mechanical  and  electrical  loading  rates.  This  is  manifested  in 
the  rate-dependent  hysteresis  behavior  observed  in  Figure  3  and  illustrates  that  hysteresis  involves 
multiple  timescales  as  well  as  spatial  scales. 

A  variety  of  techniques  have  been  employed  to  average,  or  homogenize,  grain-level  relations  to 
provide  macroscale  constitutive  relations  —  e.g.,  see  [33,51].  This  includes  Reuss  approximations, 
based  on  the  assumption  that  stress  and  electrical  fields  are  homogeneous  throughout  the  polycrys¬ 
tal  [36,60],  and  various  self-consistent  averaging  techniques  [34,35].  In  the  homogenized  energy 
framework,  this  is  accomplished  by  assuming  that  quantities  such  as  effective  and  interaction  fields 
are  manifestations  of  underlying  densities  that  are  subsequently  estimated  through  fits  to  measured 
macroscopic  data.  The  assumptions  and  techniques  used  to  homogenize  from  the  grain  to  macro-level 
often  determine  the  accuracy  and  efficiency  of  the  resulting  model. 

We  note  that  this  summarizes  only  a  few  of  the  issues  addressed  by  micromechanical  models 
and  other  research  has  focused  on  the  quantification  of  field  and  stress  interactions  between  grains 
[73],  the  incorporation  of  friction  effects  in  domain  wall  movement  [74],  micromechanics  models 
based  on  irreversible  thermodynamics  principles  [38,83],  and  direct  finite  element  implementation  of 
micromechanics  models  [37]. 

In  this  summary,  we  have  neglected  the  discussion  of  atomistic  models  that  quantify  material 
behavior  at  the  level  of  the  unit  cell.  This  comprises  a  critical  research  area  for  understanding 
fundamental  material  properties  and  for  designing  new  materials.  However,  it  is  generally  too  com¬ 
putationally  intensive  for  direct  macroscopic  material  characterization  so  we  instead  refer  the  reader 
to  [11]  for  an  overview  of  associated  models. 

Phenomenological  Models 

The  goal  when  developing  phenomenological  models  is  to  reduce  complexity  and  often  ensure  ther¬ 
modynamic  consistency  by  constructing  appropriate  relations  for  d,  Pjrr  and  Eirr  based  on  attributes 
of  measured  data.  The  reader  is  referred  to  [42,51,65]  for  overviews  of  certain  phenomenological 
frameworks. 

In  a  number  of  models,  the  remanent  or  irreversible  polarizations  and  strains  are  expressed  in 
terms  of  nonphysical  internal  variables,  often  chosen  to  satisfy  tenets  of  irreversible  thermodynamics 
[19,44]  or  emulate  theory  of  metal  plasticity  [45,52].  A  number  of  these  models  have  led  to  successful 
finite  element  implementation  [42,43,50,56,59]. 

More  recently,  investigations  have  focused  on  representing  d ,  Pirr  and  e,rr  in  terms  of  Preisach 
operators  [41],  Jiles-Atheron  hysteresis  relations  [27]  and  hysteretic  recurrent  neural  networks  [57]. 
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The  Prandtl-Ishlinskii  (PI)  hysteresis  operator  is  also  receiving  significant  attention  since  it  can  be 
proven  that  its  inverse  is  also  a  Prandtl-Ishlinskii  operator  [14] .  Whereas  the  classical  PI  operator,  like 
the  Preisach  operator  for  which  it  is  a  subset,  is  symmetric  and  rate- independent,  recent  extensions 
to  the  theory  address  asymmetry,  rate-dependence  and  certain  creep  phenomena  [2,3,48].  Although 
these  extensions  complicate  the  construction  of  inverse  models,  the  PI  model,  along  with  the  Preisach 
and  homogenized  energy  frameworks,  provides  a  feasible  avenue  for  feedforward  control  designs 
based  on  inverse  models.  Finally,  we  note  that  both  the  Preisach  and  PI  models  can  be  viewed  as 
multiscale  models  in  which  kernels,  or  hysterons,  phenomenologically  describe  grain-level  behavior 
and  homogenization  to  macroscales  is  achieved  through  integration  against  density  functions. 

Whereas  there  exist  a  number  of  rate-independent  models,  there  are  substantially  fewer  models 
that  are  capable  of  modeling  rate-dependent  phenomena  such  as  that  shown  in  Figure  3.  The 
extended  Preisach-PI  formulation  [48]  quantifies  certain  rate-dependent  material  behavior  as  does 
the  variational  formulation  [65]  and  models  of  [5,9]  which  incorporate  kinetic  theory  to  characterize 
rate-dependent  loading  effects.  The  homogenized  energy  model  (HEM),  which  is  the  topic  of  this 
paper,  quantifies  rate-dependent  effects  based  on  evolution  equations  in  combination  with  the  theory 
of  thermally  activated  processes. 

Homogenized  Energy  Model  (HEM) 

The  homogenized  energy  model  is  a  multiscale,  microscopically-motivated  or  nricronrechanical 
approach  in  the  sense  that  it  begins  with  energy  analysis  at  the  domain  level  to  construct  local 
constitutive  relations  (2).  To  construct  the  grain-level  expressions 

P(E,  a)  =  a  J2  daxa(E,  a)  +  X°E  +  P%xa(E,  a) 

a  a 

e(E,  a)  =  sEc t  +  E  ^  daxa(E ,  a)  +  ^  e%xa(E,  cr), 

a  a 

the  internal  variables  xa  are  chosen  to  be  the  fraction  of  a-variant  dipoles;  e.g.,  a  =  ±180,90  for 
tetragonal  materials.  The  dynamics  of  xa  are  governed  by  evolution  equations  driven  by  likelihoods 
constructed  using  Boltzmann  theory  to  quantify  the  scaled  probability  of  transitioning  between  stable 
equilibria  associated  with  dipole  variants.  This  incorporates  the  rate-dependence  and  multiple-time 
scales  exhibited  by  the  data  in  Figure  3. 

For  homogeneous  single  crystals  with  negligible  interaction  fields,  these  relations  can  adequately 
quantify  macroscale  behavior.  For  polycrystalline  materials  and  single  crystals  with  non-negligible 
interaction  fields,  macroscale  models  are  constructed  by  assuming  that  properties  such  as  coercive 
fields,  critical  driving  forces,  and  interaction  fields  are  manifestations  of  underlying  densities  rather 
than  constants.  This  yields  a  homogenized  energy  framework  that  accurately  characterizes  a  range 
of  material  behavior  while  retaining  the  efficiency  required  for  data-driven  parameter  estimation, 
uncertainty  quantification,  design,  and  real-time  model-based  control  implementation. 

The  framework  was  proposed  in  [75,78-80]  for  180°  ferroelectric  switching  and  extended  in  [7,71] 
to  include  90°  ferroelastic  switching.  In  both  cases,  transition  likelihoods  were  constructed  using  the 
theory  of  thermally  activated  processes  as  originally  applied  to  shape  memory  alloys  in  [1,72]  and 
magnetic  materials  in  [76,77].  It  is  shown  in  [75,79]  that  it  thus  provides  a  unified  framework  for 
characterizing  hysteresis  in  ferroic  materials.  This  approach  has  a  strong  theoretic  basis  and  provides 
accurate  and  efficient  characterization  of  1-D  ferroelectric  switching  but  becomes  increasingly  slow 
to  implement  for  ferroelastic  and  3-D  ferroelectric  switching  due  to  complications  associated  with 
computing  inflection  points  and  curves  on  2-D  and  3-D  energy  landscapes  having  4-6  local  minima.  To 
address  this,  York  used  the  driving  force  between  local  minima  to  construct  mesoscale  kernels  for  1-D 
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PZT  with  ferroelastic  switching  [92]  whereas  Kim  and  Seelecke  employed  this  technique  to  construct 
a  3-D  electromechanical  model  [47].  York  incorporated  certain  material  and  field  nonhomogeneities 
within  the  hysteron  and  Kim  modeled  certain  polycrystalline  behavior  using  a  representative  volume 
element  (RVE)  model  [46].  However,  the  combination  of  the  driving  force  likelihood  relations  with 
stochastic  homogenization  to  construct  general  macroscopic  models  for  polycrystalline  materials  has 
not  been  previously  investigated. 

The  novel  contributions  and  organization  of  the  paper  can  be  summarized  as  follows.  In  Section  2, 
we  provide  a  comprehensive  discussion  of  the  dipole  processes  that  produce  hysteresis  and  constitutive 
nonlinearities  in  ferroelectric  materials.  Whereas  much  of  this  material  is  classical,  we  focus  on 
providing  a  careful  discussion  regarding  the  sources  and  mechanisms  leading  to  even-powered  (e.g, 
quadratic)  field-strain  effects  since  this  is  important  for  model  development  and  is  often  ambiguous 
or  contradictory  in  the  literature.  In  Section  3,  we  summarize  the  180°  polarization  model  that 
was  originally  proposed  in  [78,80].  The  novel  component  of  this  section  is  Section  3.2  in  which 
we  rigorously  derive  likelihood  relations  based  on  Boltzmann  theory.  This  motivates  alternative 
likelihood  formulations  based  on  activation  energies  and  thermodynamic  driving  forces.  Section  4 
contains  a  derivation  of  the  strain-polarization  models  based  on  both  180°  and  non-180°  switching. 
The  grain-level  analysis  is  similar  to  that  in  [47,92]  but  the  development  of  macroscopic  models  based 
on  the  homogenized  energy  framework  is  new.  In  Section  4.2,  we  show  that  the  framework  yields 
constitutive  relations  of  the  form  (3)  where  d(E,a),Pirr(E,a)  and  Sirr(E,a)  are  homogenized  or 
averaged  relations  employing  the  kernel  expressions  (4).  This  framework  is  then  used  in  Section  5.1 
to  derive  a  lumped  actuator  model  as  well  as  constitutive  relations  for  distributed  systems.  The 
performance  of  the  model,  for  characterizing  creep  behavior  and  rate-dependent  phenomena  such  as 
that  shown  in  Figure  3,  is  demonstrated  in  Section  6. 

It  is  shown  in  the  companion  paper  [31]  that  due  to  its  energy  basis,  the  model  admits  highly 
efficient  implementation  and  data-driven  algorithms  to  determine  initial  parameter  estimates  based 
on  measured  properties  of  the  data.  This  facilitates  model  calibration  and  implementation  for  design 
and  control  of  devices  and  complex  structures  arising  in  applications.  The  model’s  feasibility  for 
control  applications  is  further  bolstered  by  the  fact  that  robust  inverse  model  algorithms  can  be 
implemented  at  rates  that  are  proven  no  slower  than  1/6- 1/7  the  rate  of  forward  algorithms.  Finally, 
the  unified  nature  of  the  framework  facilitates  its  extension  to  magnetic,  shape  memory  alloy,  or 
hybrid  systems. 

2  Ferroelectric  Materials 

As  detailed  in  Section  1,  commonly  employed  transducer  materials  include  BaTiC>3,  PZT,  PLZT, 
PMN  and  PVDF.  Whereas  the  specific  molecular  mechanisms  are  material-dependent,  all  of  these 
compounds  exhibit  certain  shared  rneso-  and  macro-scale  properties  which  form  the  basis  for  homog¬ 
enized  energy  models.  In  this  section,  we  summarize  relevant  material  properties;  we  refer  the  reader 
to  [40, 58, 66,  75]  for  details. 

Barium  titanate  and  lead  titanate  (PbTiOs)  are  isostructural  with  the  mineral  perovskite  (CaTiOs) 
and  exhibit  what  is  termed  a  ABO3  perovskite  structure.  As  illustrated  in  Figure  5  for  BaTiC>3,  this 
consists  of  a  paraelectric  non-polar  cubic  structure  above  the  Curie  point  Tc  and  ferroelectric  tetrag¬ 
onal,  orthorhombic,  and  rhombohedral  forms  at  temperature  below  Tc.  It  is  noted  that  for  typical 
operating  temperatures,  BaTiC>3  and  PbTiC>3  exhibit  tetragonal  structures  whereas  the  composition 
8/65/35  PLZT  (Pbo.08Ti.35Zr.65O3)  is  near  a  rhombohedral-tetragonal  morphotropic  phase  boundary 
and  PbZr03  is  orthorhombic.  We  will  focus  primarily  on  tetragonal  structures  but  analogous  results 
hold  for  rhombohedral  and  orthorhombic  forms. 
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Cubic  Tetragonal  Orthorhombic  Rhombohedral 


Figure  5:  Cubic,  tetragonal,  orthorhombic  and  rhombohedral  forms  of  perovskite  compounds  and 
approximate  transition  temperatures  for  BaTiC>3  (°C). 

As  illustrated  in  Figure  6  for  BaTiOs,  the  materials  exhibit  a  spontaneous  polarization  Po  f°r 
T  <  Tc.  This  leads  to  to  180°  ferroelectric  switching  for  fields  greater  than  a  coercive  field  Ec  and 
90°  ferroelastic  switching  for  compressive  stresses  larger  than  a  coercive  stress  ac.  Both  effects  can 
occur  in  transducers  and  for  devices,  such  as  THUNDER,  that  exhibit  large  internal  prestresses, 
both  effects  may  be  significant. 

The  electromechanical  properties  of  ferroelectric  compounds  are  intimately  related  to  the  manner 
in  which  the  polar  tetragonal,  orthorhombic  or  rhombohedral  structures  respond  to  input  fields  and 
stresses.  To  illustrate  these  effects,  we  summarize  the  polarization  and  strain  behavior  for  three  cases: 
(i)  single  crystal,  180°  ferroelectric  domains,  (ii)  single  crystal,  180°  ferroelectric  and  90°  ferroelastic 
domains,  and  (iii)  polycrystalline  materials. 

Single  Crystal,  180°  Ferroelectric  Domains 

In  the  absence  of  applied  or  internal  stresses,  minimization  of  the  electrostatic  energy  yields 
twinned  180°  domains.  For  fields  E  less  in  magnitude  than  the  coercive  field  Ec,  this  yields  ap¬ 
proximately  linear  E-P ,  E-e,  P-e ,  where  e  is  the  strain,  behavior  as  illustrated  in  Figure  7.  This 
is  due  to  reversible  deformations  of  the  tetragonal  cell  and  resulting  strains  are  small  compared  to 
those  resulting  from  90°  switching.  Fields  in  excess  of  Ec  produce  180°  switching  which  yields  large 
changes  in  polarization  but  small  changes  in  strains  since  differences  in  the  latter  are  due  only  to 
small  changes  in  the  cell  dimension. 


(a)  (b)  (c) 


Figure  6:  (a)  Tetragonal  structure  of  PbTiOs  for  T  <TC  and  resulting  spontaneous  polarization  Po- 
(b)  Ferroelectric  180°  polarization  switch  due  to  an  applied  electric  field  E  >  Ec  and  (c)  ferroelastic 
90°  degree  switch  due  to  a  compressive  force  a  larger  in  magnitude  than  the  coercive  stress  ac. 
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Figure  7:  (a)  Polarization  changes  due  to  180°  ferroelectric  dipole  switching,  and  (b),  (c)  small 
strains  resulting  from  deformations  of  the  unit  cell. 

Single  Crystal,  180°  Ferroelectric  and  90°  Ferroelastic  Domains 

The  situation  is  more  complicated  when  fields  and  stresses  are  simultaneously  applied  to  the 
materials.  We  illustrate  the  case  of  compressive  stresses  since  this  induces  dipole  switching.  However, 
transducers  such  as  macro-fiber  composites  (MFC)  and  THUNDER  will  exhibit  regions  having  both 
tensile  and  compressive  stresses. 

We  first  note  that  the  strains  resulting  from  90°  switching  are  significantly  larger  than  the  linear 
strains  resulting  from  reversible  deformations  of  the  unit  cell.  For  BaTiC>3  single  crystals,  this  90° 
switching  can  produce  strains  up  to  1.1%.  The  reader  is  referred  to  [15]  for  experimental  results 
illustrating  the  E-P,  E-e  and  P-e  behavior  of  single  crystal  BaTiC>3  subject  to  various  prestress 
levels.  The  behavior  for  a  =  —1.78  MPa  is  plotted  in  Figure  ll(a)-(c). 

In  Figure  8,  we  illustrate  ferroelastic  switching  for  an  ideal  single  crystal.  As  compressive  stresses 
are  increased,  two  mechanisms  produce  changes  in  the  strains  and  polarization:  (i)  linear  elastic 
changes  in  the  tetragonal  cell,  and  (ii)  nucleation  and  growth  of  90°  domains.  The  primary  90° 
rotation  occurs  at  the  coercive  stress  ac  after  which,  the  elastic  stiffness  c  =  is  approximately 
linear.  Representative  domain  configurations  are  shown  in  Figure  8(c). 


(c) 


Figure  8:  (a)  and  (b)  Elastic  response  and  stress-induced  ferroelastic  switch  at  ac  for  E  =  0. 
(c)  Nucleation  and  growth  of  90°  domains  during  the  switching  process. 
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The  E-P,  E-e  and  P-e  behavior  at  the  points  B  and  F  corresponding  to  fixed  stresses  a  £  (— ac,  0) 
and  a  £  (—00,  —  ac)  are  depicted  in  Figure  9.  A  comparison  between  Figure  9(a)  and  Figure  7  (a  =  0) 
illustrates  that  for  a  £  (— <tc,0),  the  magnitude  of  the  polarization  remains  approximately  the  same 
but  the  behavior  near  coercivity  reflects  the  90°  switching.  The  strains  are  significantly  larger  due 
to  the  90°  switching  as  demonstrated  by  experimental  BaTiC>3  data  in  [15].  This  regime  provides 
optimal  performance  for  single  crystal  actuators. 

For  fixed  prestresses  larger  in  magnitude  than  ac,  ferroelastic  switching  yields  90°  domains  as  a 
starting  state  for  E  =  0  and  the  material  behaves  in  a  linear  elastic  manner  until  field-induced  90° 
switching  occurs  at  the  point  ii  shown  in  Figure  9(b).  Due  to  the  magnitude  of  the  prestresses,  large 
fields  E  are  required  to  produce  significant  polarization  and  strain  levels  and  single  crystal  actuators 
operating  in  these  regimes  would  have  diminished  outputs  —  or  would  be  destroyed  by  internal  field 
and  stress  levels. 

Polycrystalline  Materials 

Whereas  single  crystal  BaTiC>3  is  being  considered  for  applications,  most  ferroelectric  transducer 
materials,  such  as  PZT,  are  polycrystalline  which  affects  their  behavior  in  a  variety  of  ways.  It  is 
first  noted  that  whereas  a  single  crystal  is  polar  for  T  <  Tc,  polycrystalline  compounds  exhibit  zero 


(b) 

Figure  9:  E-P,  E-e  and  P-e  behavior  at  fixed  stresses  a  corresponding  to  the  points  (a)  B  and  (b)  F 
from  Figure  8.  Field  levels  in  (b)  are  larger  than  those  in  (a). 
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(d) 


Figure  10:  (a)  Unpoled  material,  (b)  rotation  of  180°  domains  due  to  poling,  and  (c)  90°  switching 
due  to  an  applied  field,  (d)  Ferroelastic  90°  switching  due  to  an  applied  stress. 


net  polarization  due  to  the  random  orientation  of  grains  and  domains.  This  necessitates  that  the 
materials  be  poled  before  use  as  illustrated  in  Figure  10(b)  and  (c).  Due  to  the  random  orientation 
of  grains,  the  materials  exhibit  smaller  remanence  polarizations  than  corresponding  single  crystals. 
Moreover,  as  shown  in  Figure  10(c),  polycrystalline  materials  exhibit  field-induced  90°  switching  in 
the  absence  of  a  prestress  and  both  90°  and  180°  switching  can  occur  at  stresses  and  fields  below  single 
crystal  coercive  values.  This  causes  polycrystalline  E-P ,  E-e  and  P-e  hysteresis  curves  to  be  smoother 
than  their  single  crystal  counterparts  and  it  motivates  the  use  of  models  such  as  the  homogenized 
energy  model,  Preisach  formulations  or  Prandtl-Ishlinskii  models  that  can  accommodate  the  random 
grain  orientations  via  densities  in  the  representation. 

Poly  crystalline  PLZT  data,  corresponding  to  the  ideal  single  crystal  behavior  depicted  in  Fig¬ 
ures  8  and  9,  can  be  found  in  [61].  Single  crystal  BaTiC>3  data  from  [15],  polycrystalline  PZT  data 
from  [92],  PLZT  data  from  [69],  and  PMN  data  from  [29]  are  plotted  in  Figure  11.  It  is  first  noted 
that  polycrystalline  switching  behavior  is  smoother  than  for  BaTi03  single  crystal  due  to  material 
nonhomogeneities.  It  is  also  observed  that  PLZT  exhibits  nearly  quadratic  D-E  behavior  and  PMN 
exhibits  nearly  quadratic  and  anhysteretic  E-e  behavior  as  compared  with  the  butterfly  curves  mea¬ 
sured  for  PZT.  We  provide  further  discussion  regarding  the  piezoelectric,  quadratic,  and  switching 
behavior  of  these  compounds  in  the  next  subsection. 

2.1  Piezoelectric,  Electrostrictive  and  Domain  Switching  Behavior 

Ferroelectric  materials  exhibit  a  complex  combination  of  linear,  quadratic  or  quartic,  and  hysteretic 
behavior  in  general  operating  regimes.  These  effects  are  often  categorized  as  piezoelectric,  elec¬ 
trostrictive,  or  domain  switching  in  nature  but  there  is  significant  interplay  between  the  underlying 
mechanisms  and  hence  ambiguity  in  the  definitions  of  the  phenomena.  We  detail  aspects  of  these 
phenomena  to  clarify  the  underlying  physical  mechanisms  and  motivate  issue  that  must  be  addressed 
in  models. 

The  direct  piezoelectric  effect  constitutes  the  change  in  polarity  that  results  from  an  applied 
stress  whereas  the  converse  piezoelectric  effect  constitutes  reversible  strains  generated  by  applied 
field.  In  both  cases,  the  designation  piezoelectric  has  a  linear  connotation. 
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Figure  11:  (a)  E-P,  (b)  E-e  and  (c)  P-e  behavior  of  prestress( 
(d)  E-P ,  (e)  E-e,  and  (f)  P-e  behavior  of  PZT  data  from  [92], 
of  PLZT  data  from  [69].  (j)  E-P  and  (k)  E-e  behavior  of  PIN 
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lgle  crystal  BaTi03  data  from  [15]. 
E-D,  (h)  E-e,  and  (i)  D-e  behavior 
ata  from  [29]. 


Electrostriction  is  classically  used  to  designate  field-induced  strains  that  are  proportional  to 
even  powers  of  the  field  and  hence  are  independent  of  its  polarity.  Classical  electrostriction  results 
from  reversible,  and  hence  anhysteretic,  deformations  in  the  lattice  structure  and  it  is  inherent 
to  all  materials  including  gaseous,  liquid,  and  crystalline  or  amorphic,  polar  or  centrosymmetric 
solids  [16,40,58,66].  In  most  ferroelectric  materials,  electrostrictive  effects  are  small  compared  to 
linear  piezoelectric  effects  and  are  hence  neglected.  However,  in  some  compounds  such  as  PMN,  it  is 
significant  and  must  be  included  in  models.  Depending  on  the  independent  variables,  electrostriction 
can  be  modeled  as  a  quadratic,  or  higher-order,  relation 

e  =  ME2  ,  e  =  QP2  ,  £  =  QD2 

between  fields  E,  polarization  P,  or  D-fields  and  the  resulting  strain  e. 

A  source  of  confusion  stems  from  the  fact  that  the  quadratic  dependence  of  strains  on  P,EoiD  is 
due  to  two  different  mechanisms:  classical  electrostriction  and  domain  rotation.  As  noted  previously, 
electrostriction  results  from  reversible  deformations  of  the  lattice  structure;  hence  it  is  anhysteretic 
and  it  occurs  to  varying  degrees  in  all  materials.  In  ferroelectric  materials,  domain  rotation  only 
occurs  for  T  <TC  when  the  material  is  in  a  polar  state.  The  quadratic  strains  due  to  domain  rotation 
are  typically  larger  than  electrostrictive  effects,  except  in  materials  such  as  PMN,  and  E  —  e  curves 
exhibit  hysteresis.  The  difference  between  the  nearly  anhysteretic,  electrostrictive,  behavior  of  PMN 
and  the  hysteretic,  domain  switching,  behavior  of  PZT  and  PLZT  is  illustrated  in  Figure  11. 

We  note  that  some  researchers  define  quadratic  effects  due  to  domain  reorientation  as  electrostric¬ 
tion  [84]  whereas  others  more  generally  define  electrostriction  as  field-induced  strains  that  are  in¬ 
dependent  of  the  field  polarity.  We  follow  instead  the  distinction  made  by  Caspari  and  Merz  [17] 
between  classical  electrostriction  and  quadratic  domain  rotation  effects  since  the  two  are  character¬ 
ized  differently  in  the  homogenized  energy  model. 

We  next  detail  why  domain  reorientation  produces  quadratic  strain  dependencies  and  then  illus¬ 
trate  the  relative  contributions  of  electrostriction  and  quadratic  domain  reorientation  for  common 
ferroelectric  materials. 

Quadratic  Strains  Due  to  Domain  Rotation 

To  illustrate  how  domain  rotation  can  produce  quadratic  strains,  we  represent  paraelectric  cubic 
regions  by  unit  spheres  and  the  ferroelectric  polar  regions  by  ellipsoids  as  motivated  by  analysis  of 
magnetostriction  presented  on  pages  343-346  of  [18]  or  pages  99-104  of  [39].  We  let  e  denote  the 
spontaneous  strain  generated  by  the  paraelectric  to  ferroelectric  phase  transition  when  all  dipoles 
are  aligned  with  the  field.  As  depicted  in  Figure  12(a),  the  strain  at  the  angle  8  is  thus 

s(8)  =  e  cos2  9.  (5) 

For  a  random  domain  orientation,  the  spontaneous  strain  for  the  polarized  material  is  thus 

Z*71'/2  g: 

Ao  =  /  £  cos2  6  sin  9d9  =  - . 

Jo  3 

The  application  of  a  field  causes  domains  to  rotate  in  the  manner  depicted  in  Figure  12(c)  thus 
yielding  the  relative  change  in  dimension  As  =  £  —  Ao  =  |e.  If  we  assume  that  changes  in  polarization 
are  solely  due  to  rotation,  substitution  of  the  relation  P  =  PqcosO  and  e  =  §Ao  into  (5)  yields  the 
quadratic  strain  relation 

e(P)  =  § 
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Figure  12:  (a)  Spontaneous  strain  as  a  function  of  the  angle  6.  (b)  Spheres  used  to  model  isotropic, 
disordered  material  behavior  in  the  paraelectric  phase  and  (c)  ellipsoids  representing  the  order  en¬ 
capsulated  by  domains  for  T  <  Tc.  (d)  Total  strain  e  and  strain  As  due  to  orientation  of  domains  in 
the  field  direction. 


2.2  Material  Behavior 
Magnesium  Niobate  Compounds 

As  detailed  in  [75],  PMN-PT-BT  transducers  are  employed  in  applications  including  sonar  trans¬ 
duction  due  to  their  low  hysteresis,  high  strain  capabilities.  These  relaxor  ferroelectric  compounds 
also  exhibit  significant  electrostrictive  behavior  as  illustrated  in  Figure  1 1  (j )- (k)  with  data  from  [29]. 
For  an  operating  temperature  of  5°  C,  the  material  exhibits  nearly  quadratic  E-e  behavior  at  low  to 
moderate  field  inputs  and  even,  higher-order  behavior  at  high  fields.  This  provides  a  benchmark  for 
discussing  the  degree  to  which  electrostriction  plays  a  role  in  BaTiC>3,  PZT  and  PLZT  behavior. 

Single  Crystal  Barium  Titanate 

The  behavior  of  BaTiOs  has  been  heavily  studied  for  over  60  years  and  it  remains  an  important 
research  topic  due  to  the  capability  of  single  crystals  to  generate  large  strains  (~1  %)  through 
90°  rotations.  The  behavior  of  single  crystal  BaTiOs  subjected  to  various  compressive  stresses  has 
recently  been  investigated  in  [15]  and  representative  results  are  shown  in  Figure  ll(a)-(c).  It  is 
observed  that  in  contrast  to  PMN-PT-BT,  the  BaTiOs  E-e  curve  exhibits  significant  hysteresis 
which  gives  it  a  classical  butterfly  profile.  Early  authors  attributed  this  to  electrostrictive  effects  [63] 
but  analysis  originating  in  the  classic  1950  paper  by  Caspari  and  Merz  [17]  demonstrates  that  it  is 
primarily  due  to  domain  reorientation.  This  is  the  source  of  the  hysteresis. 

In  [17],  Caspari  and  Merz  show  that  there  are  three  sources  of  quadratic  strain  effects  in  BaTiOs: 
spontaneous  strains  e o  =  QPq  resulting  from  the  paraelectric  to  ferroelectric  phase  transition,  but¬ 
terfly  effects  due  to  dipole  switching,  and  classical  electrostrictive  effects.  To  analyze  the  relative 
contribution  due  to  electrostriction,  they  consider  the  field-induced  strains  e  that  result  when  the 
polarization  is  biased  about  the  spontaneous  polarization  P$  and  corresponding  field  Eq  as  depicted 
in  Figure  13.  Application  of  an  external  field  yields  a  change  As  in  the  strain  so  that 

e  =  £q  +  As  =  Q{Pq  +  P )2 

=>  As  =  {2QP0)P  +  QP2. 

If  we  take  P  =  yE,  the  change  in  strain  can  be  expressed  as 

As  =  (2QXP0)E  +  ( QX2)E 2. 

It  follows  that  the  piezoelectric  constants  g  and  d  and  the  electrostrictive  constant  M  are  given  by 

g  =  2 QPq  ,  d  =  2QXPo  ,  M  =  Qy2. 
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Figure  13:  Change  in  strain  Ae  due  to  applied  field  increments  about  Eq  =  Pq/x- 


For  the  cgs  parameter  values  given  in  [17],  it  is  unfortunately  difficult  to  directly  compute  the 
relative  strain  contributions.  Instead,  we  use  the  following  parameter  values  from  [8]  —  see  also  [55]  - 

Q  =  0.11  m4/C2 

P0  =  0.26  C/m2 

X  =  eoXe  =  9.56  x  1CT10  C/(Vm)  ,  Xe  =  er  -  1  ,  er  =  109 

which  yields 

g  =  5.72  x  10~2  m2/C  ,  [8]  report  5.75  x  10-2  m2/C  measured 

d  =  5.47  x  1CT11  m/V  ,  [8]  report  8.56  x  1CP11  m/V  measured 

M  =  1.01  x  1CP19  m2/V2. 

For  a  change  in  polarization  P  =  0.05  C/m2  and  field  E  =  5  x  105  V/m,  which  is  representative  of 
the  data  in  [15],  this  yields  the  strain  contributions 

gP  =  2.86  x  1CT3  ,  QP2  =  2.75  x  10~4 
dE  =  2.74  x  10~5  ,  ME 2  =  2.53  x  10~8. 

In  both  cases,  the  quadratic  electrostrictive  effects  are  negligible  compared  with  the  linear  piezo¬ 
electric  effects  and  hence  can  be  neglected.  This  implies  that  the  piezoelectric  effect  can  be  interpreted 
as  the  electrostrictive  effect  biased  about  the  spontaneous  polarization  Pq.  It  further  implies  that 
the  E-e  and  P-e  behavior  shown  in  Figure  ll(a)-(c),  which  is  independent  of  field  polarity,  is  due 
to  the  90°  and  180°  switching  which  imbues  the  material  with  large  strain  capabilities.  The  P-e 
behavior  in  Figure  11(c)  also  demonstrates  that  180°  and  90°  switching  for  BaTiC>3  occur  on  very 
different  timescales.  The  flat  region  indicates  that  180°  polarization  switching  occurs  in  advance 
of  strain-producing  90°  switches.  It  is  shown  in  Section  3.3  and  [31]  that  single  crystal  behavior 
of  this  type  can  be  modeled  using  energy  relations  that  incorporate  both  90°  and  180°  switching 
mechanisms  and  employ  linear  domain-level  constitutive  relations. 

Polycrystalline  Tetragonal  PZT  and  PL8ZT 

As  illustrated  in  Figure  10,  polycrystalline  compounds  exhibit  90°  switching  in  the  absence  of 
applied  stresses  and  at  low  field  levels.  To  quantify  the  role  of  90°  switching  in  poly  crystalline 
compounds,  Tsurumi  et  al.  employed  X-ray  diffraction  (XRD)  techniques  to  measure  the  amount 
of  90°  reorientation  [84].  We  note  that  the  PL8ZT  compounds  under  consideration  had  chemical 
compositions  which  ensured  that  they  were  in  the  tetragonal  form  and  not  near  the  morphotropic- 
phase  boundary  (MPB).  This  is  in  contrast  to  8/65/35  PLZT  which  is  at  a  tetragonal-rhombohedral 
nrorphotropic  boundary. 
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Figure  14:  Measured  strains  emeas,  strains  ep  =  d^E  due  to  linear  piezoelectric  effects,  and  strains 
ep+e 90  due  to  piezoelectric  effects  and  90°  dipole  rotation  for  (a)  PZT  and  (b)  PL8ZT  data  from  [84], 


Biased  minor  strain  loops  for  PZT  and  PL8ZT  data  from  [84]  are  shown  in  Figure  14  and  the 
P-E  hysteresis  curve  and  E-e  butterfly  curve  for  PL8ZT  are  plotted  in  Figure  15.  In  each  case,  emeas 
is  the  measured  strain,  ep  =  d^E  is  the  computed  piezoelectric  strain,  and  ego  is  the  strain  due  to 
90°  dipole  rotation  as  measured  using  the  XRD  techniques.  Figure  14  illustrates  that  piezoelectric 
strains  contribute  less  than  half  of  the  measured  response  and  that  the  sum  of  the  piezoelectric 
and  rotational  effects  are  consistent  with  the  measured  strains.  The  significant  contribution  due 
to  90°  reorientation  is  further  demonstrated  for  PL8ZT  in  Figure  15.  Here  it  is  observed  that  90° 
reorientation  constitutes  a  primary  strain  contribution  at  fields  below  the  coercive  field. 

These  results  demonstrate  that  90°  reorientation  is  one  of  the  primary  strain  producing  mech¬ 
anisms  in  polycrystalline  materials,  especially  at  low  field  levels.  Hence  it  must  be  incorporated  in 
models  to  achieve  accurate  material  characterization. 

Furthermore,  Li  et  al.  [54]  demonstrate  using  XRD  analysis  that  polarization  switching  in  PZT 
compounds,  with  compositions  near  the  morphotropic  boundary,  are  due  primarily  to  two  successive 
90°  switches  rather  than  a  single  180°  transition.  It  is  further  noted  in  [96]  that  180°  switching 
occurs  more  rapidly  than  90°  switching  as  exemplified  for  BaTiC>3  by  the  nearly  flat  P-e  region  in 
Figure  11(c).  Hence  we  model  180°  switches  with  a  different  timescale  than  the  90°  rotations  that 
produce  large  changes  in  strains. 


(a)  (b) 

Figure  15:  (a)  E-P  PL8ZT  data  from  [84],  (b)  Total  field- induced  strains  £meas  and  strains  ego  due 
to  90°  dipole  rotation  for  PL8ZT  data  from  [84], 
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Rhombohedral  PLZT 


The  PLZT  compound  8/65/35  differs  from  the  previously  described  PL8ZT  material  in  that  it 
has  a  rhombohedral  rather  than  tetragonal  crystal  structure.  Hence  polarization  switches  of  70.5° 
(three  nearest  corners)  or  109.5°  (three  further  corners)  produce  deformations  or  strains  in  the 
material  [9,61].  Like  PMN,  it  is  also  a  relaxor  ferroelectric  material. 

It  is  illustrated  in  Figure  11  (g)-(i)  that  the  D-e  behavior  of  8/65/35  PLZT  is  quadratic  and  nearly 
anhysteretic  whereas  the  E-e  curve  has  the  usual  hysteretic  butterfly  shape.  In  [94],  these  effects 
were  modeled  as  electrostrictive  strains.  However,  they  are  more  likely  due  to  domain  switching  for 
the  previously  mentioned  reasons. 

There  are  various  options  to  model  the  quadratic  strain  behavior.  One  is  to  construct  energy 
functionals  having  quadratic  polarization  or  D-field  dependence  [69] .  This  models  the  quadratic  P-e 
or  D-e  behavior  in  a  phenomenological  sense.  A  second  alternative  is  to  employ  energy  relations 
that  incorporate  the  dipole  switching  mechanisms  in  combination  with  linear  constitutive  relations. 
This  is  the  strategy  that  we  employ  in  Section  3  and  4. 


3  Homogenized  Energy  Model  —  Polarization 

The  homogenized  energy  model  (HEM)  for  ferroelectric  materials  is  a  multiscale  approach  comprised 
of  two  fundamental  components:  (i)  construction  of  energy-based  grain-level  kernels  that  characterize 
dipole  switching,  material  constitutive  behavior,  and  thermal  relaxation  mechanisms,  and  (ii)  con¬ 
struction  of  macroscale  models  through  the  assumption  that  coercive  fields,  critical  driving  forces, 
and  interaction  fields  are  manifestations  of  underlying  densities  rather  than  constants. 

Ferroelectric  materials  subjected  to  general  fields  E  =  [E\ ,  ±2,  E:$]  or  stresses  a  exhibit  vector¬ 
valued  polarization  P  =  [P\,  P2,  Ps\  polarization  or  strain  e  responses.  As  illustrated  in  Section  1, 
however,  many  actuators  and  sensors  employ  1-D  input  and  output  responses  so,  for  these  applica¬ 
tions,  it  is  advantageous  to  consider  scalar  electric  and  mechanical  variables  E  =  E$,  a  =  <733,  P  =  P3 
and  e  =  £33.  For  applications  that  require  truly  2-D  or  3-D  polarization  or  strain  relations,  energy 
functionals  can  be  constructed  using  the  framework  in  [46, 47] . 

We  summarize  first  the  homogenized  energy  model  for  polarization.  As  illustrated  in  Figure  7, 
the  primary  polarization  features  for  small  force  regimes  can  be  defined  by  considering  only  ±180° 
dipole  states.  Whereas  this  case  has  been  previously  reported  in  [75,78-80],  the  summary  provided 
here  motivates  and  provides  the  context  for  the  strain-polarization  model  in  Section  4.  Moreover, 
we  present  new  theory  in  Section  3.2  for  the  transition  likelihoods. 

3.1  Polarization  Kernels  with  Negligible  Thermal  Activation 

To  characterize  180°  dipole  switching,  we  employ  the  Gibbs  energy  density1 

G_  (E,  p )  =  \p{P  +  PR )2  —  EP  ,  P< -Pi 

G+(E,P)  =  \p{P-PR)2-EP  ,P>Pi  (6) 

GU(E,  P)  =  ln(Pr  -  PR)  (g  -  PR )  -  EP  ,  |P|  <  Pj 

lrThe  fact  that  G  is  an  energy  density  results  from  the  polarization  definition  P  =  linidy^o  |  -^y  pij  where  Pi 

is  a  collection  of  general  dipoles  in  a  nonhomogeneous  region,  dV  is  a  reference  volume  and  Nv  denotes  the  number  of 
dipoles  per  unit  volume  [6] .  The  polarization  thus  designates  the  dipole  moment  density  of  the  material  and  has  units 
of  coulombs  per  square  meter  (C/m2).  This  motivates  some  authors  to  use  the  terminology  “polarization  density”  or 
“dipole  moment  density”  when  introducing  the  polarization  [23,26]. 


G(E,P)  =  ip(P)  —  EP  =  { 
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where  ip  is  a  piecewise  quadratic  Helmholtz  energy  and  Pr,  t]  and  Pj  respectively  denote  the  rerna- 
nence  polarization,  inverse  susceptibility  after  switching,  and  the  positive  inflection  point.2 

For  regimes  in  which  thermal  excitation  is  negligible,  the  condition  jjp  =  0,  which  reflects  the 
reorientation  of  dipoles  to  achieve  energy  minimization,  yields  the  linear  local  polarization  relations 

P~(E)  =  --PR,  P+(E)  =  -  +  Pr,  Pm(E )  =  PlE  (7) 

v  n  -  pR) 

The  Gibbs  energy  density  corresponding  to  these  two  equilibrium  conditions  is 

Ga{E)  = -l-XE2  -  EP%  (8) 

where  a  =  ±180  and  y  = 

/V  fj 

To  specify  a  hysteresis  kernel,  or  hysteron,  we  consider  an  ideal  lattice  with  volume  V  having  N 
cells  of  the  form  depicted  in  Figure  6(a).  We  let  JV_  and  N+  denote  the  number  of  negatively  and 
positively  oriented  dipoles  so  the  corresponding  dipole  fractions  are  x+  =  ,  X-  =  Because 

N+  ±  iV_  =  1,  it  follows  immediately  that 

x+  ±  X-  =  1.  (9) 

We  let  P~,  P+  denote  the  polarizations  due  to  negative  and  positive  dipoles  and  let  pafj,  a,  (3  =  ±180, 
which  has  units  of  (1/s),  denote  the  likelihood  of  transitioning  form  an  a- well  to  /3-well. 

The  evolution  of  dipole  fractions  is  governed  by  the  differential  equation 

x+  =  -p+-x+  +  p-+x- 

X-  =  p- 1 _ X+  —  p _ \-X- 

which  can  be  simplified  to 

X+ = -P+-X+ +p-+{  1-X+)  (10) 

using  the  identity  (9).  The  hysteron  is  given  by 


P  =  X+P++X-P  .  (11) 

The  specification  of  p _ \-,P-\ _ ,P+,P~  depends  on  the  degree  to  which  kinetics  due  to  thermal 

activation  are  relevant.  For  operating  conditions  in  which  relaxation  due  to  thermal  activation  is 
negligible,  it  is  shown  in  [75]  that  transition  likelihood  rates  are  given  by 

l  ,  E>  Ec  f  l  ,  E  <  -Ec 

0  ,  else  ’  P+~  ~  {  0  ,  else 

where  1/r  is  the  frequency  at  which  dipoles  attempt  to  switch.  As  illustrated  in  Figure  7,  Ec  denotes 
the  coercive  field  which  is  related  to  Pi  and  Pr  via  the  expression  Ec  =  i](Pr  —  Pi).  For  this  operating 
regime,  p+  =  p+  and  P  =  P,m  so  the  hysteron  is 

P  = - Pr  ±  2x.| -Pr.  (13) 

V 

To  incorporate  the  kinetics  due  to  thermal  activation,  which  produces  creep  and  accommodation 
or  reptation-like  effects,  likelihoods  can  be  specified  in  terms  of  two  quantities:  error  functions 
resulting  from  the  theory  of  thermally  activated  processes  or  activation  energies  A Ga- 

2We  employ  (6)  to  collectively  specify  the  Gibbs  energy  density  for  both  the  positive  and  negative  polarization 
variants.  Because  the  global  expression  (6)  is  bistable,  it  is  not  a  Gibbs  energy  density  in  the  classical  sense  which 
has  prompted  some  authors  to  designate  it  a  Landau  energy  [10,25].  Because,  we  are  not  employing  higher-order 
polynomials,  as  is  typically  the  case  for  the  Landau  energy,  we  follow  the  precedent  set  by  authors  such  as  Devonshire  [22] 
who  use  this  global  representation  with  the  understanding  that  its  reversible  Legendre  transform  properties  hold  for 
the  individual  variants. 
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3.2  Theory  of  Thermally  Activated  Processes 

We  illustrate  the  theory  in  the  context  of  determining  the  likelihood  that  a  dipole  switches  from 
negative  to  positive.  If  we  sum  over  the  index  set  of  negative  dipoles,  the  resulting  polarization  is 

p-  =  vT.N-p- 

i—  1 

However,  for  the  purpose  of  energy  formulation  and  minimization,  it  is  advantageous  to  sum  over 
the  set  of  possible  states  rather  than  the  index  set  For  a  finite  set  of  negative  dipole  states  S- ,  this 
yields  the  relations 

a*) 

p£S-  p£S- 

where  Np  is  a  distribution  quantifying  the  number  of  dipoles  having  strength  p  as  illustrated  in 
Figure  16(b).  For  a  continuum  of  dipole  values,  IV_  and  P~  are  given  by 

N-  =  f  N(p)dp  ,  P~  =  I  pN(p)dp.  (15) 

JpeS-  V  JPeS- 

We  consider  first  the  case  when  S-  is  finite. 

The  Gibbs  energy  for  negatively  oriented  dipoles  is 

g  =  &  +  K- ST-  EV~ 

where  $  and  K  denote  the  internal  and  kinetic  energies,  S  is  the  entropy,  T  is  temperature  in  degrees 
Kelvin,  and  V~  is  the  total  dipole  strength.  We  note  that  g  has  units  of  V-C  and  hence  it  is  an 
energy  rather  than  an  energy  density  like  G  given  in  (6).  The  internal  energy  is  taken  to  be 

*  =  J2  ^np 

p£S- 


where  cj>(p)  quantifies  the  internal  energy  for  each  dipole  of  strength  p.  Letting  po  denote  the  spon¬ 
taneous  strength  of  a  negatively  oriented  dipole,  an  appropriate  choice  for  < j>(p)  is 


<j)(p) 


1 

2 


v{p  +  po)2  +  m- 


(b) 


Figure  16:  (a)  Dipole  configurations  used  to  construct  p _ using  the  theory  of  thermally  activated 

processes  (ooo»),  the  activation  energy  AG“  (+,  •)  and  the  thermodynamic  driving  force  F _ (_  (+,  x). 

(b)  Negative  dipole  states  p  and  density  Np. 
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Following  [72],  where  analogous  arguments  are  presented  for  SMA,  we  assume  that  K  is  a  linear 
function  of  temperature  and  take 

K  =  A(T  -  Tr) 

where  Tr  is  the  reference  temperature  where  the  internal  energy  satisfies  Ur  =  0.  This  is  analogous 
to  the  molecular  kinetic  energy  relation  K  =  \kT. 

To  compute  the  entropy,  we  note  that 


IlpeS-  Np- 

quantifies  the  number  of  ways  to  arrange  1V_  dipoles  so  that  their  states  p  have  the  distribution  Np. 
Boltzmann’s  law, 

5  =  klnW, 

then  quantifies  the  entropy  which  can  be  approximated  by 

S  =  kN_  In  iV_  —  k  ^2  np  ln  np 
p£S- 

using  Stirling’s  formula  ln  x\  ~  x  ln  x  —  x.  Finally,  we  let 

v-=YJ  pnp 

p£S- 


denote  the  total  dipole  strength.  Note  that  V~  differs  from  the  polarization  P~  in  the  sense  that 
the  latter  reflects  the  dipole  strength  per  unit  volume  which  renders  it  a  density;  compare  with  (7) 
or  (15). 

The  kinetics  of  dipole  motion  are  governed  by  a  partial  equilibrium  of  the  internal,  kinetic, 
entropic  and  electrostatic  energies  subject  to  the  constraint  N-  =  J2peg_  Np.  To  reformulate  this  as 
an  unconstrained  optimization  problem,  we  employed  the  augmented  Gibbs  functional 

SS  =  9  +  ^(n_  -  ^2  Np 
^  PeS- 

=  [^(P)  +  kTlnNP  ~  EP\  np  +  A(T  ~  tr)  ~  kTN-  In  1V_  +  j( N_  -  J2  Np 

peS-  ^  peS- 


where  7  is  a  Lagrange  multiplier.  For  a  fixed  set  of  dipole  states,  the  distribution  Np  is  determined 
by  the  equilibrium  condition  5g^  =  0  where  5g 7  is  the  first  variation  with  respect  to  Np.  Under  the 
assumption  that  for  a  given  state  p,  Np  is  sufficiently  large  to  permit  consideration  of  the  derivative 
this  yields 


Np  =  exp 


1  exp 


f  0(p)  ~  EpY 

V  kT  )_■ 


From  (14),  it  follows  that 


Np  exP 

'  4>{p)—Ep\ 

<  kT  )\ 

N 

z2pgs.  exP 

(  Np)-Ep\ 

Y  V  kT  )\ 

(16) 


We  note  that  (16)  quantifies  the  probability  of  observing  Np  dipoles  for  each  dipole  state  p  as 
illustrated  in  Figure  16(b).  It  can  also  be  interpreted  as  the  probability  that  a  dipole  p,  subjected 
to  a  field  E,  has  an  energy  level 

g{E,p)=4>{p)-Ep\  (17) 
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that  is 


(18) 


The  likelihood  rate  p _ |_ 


Ms) 


exp  [-g(E,p)/kT\ 

52peS-exP[-9(E,p)/kT}' 


that  a  dipole  has  the  energy  to  switch  from  negative  to  positive  is  thus 


1  exp[-g(E,pM)/kT\ 

T  £PeS_  exP  [~9(E,p)/kT\ 


where  pm  is  the  dipole  strength  of  the  unstable  equilibrium  as  shown  in  Figure  16(b).  Note  that  1/r 
is  the  frequency  at  which  dipoles  attempt  to  switch  so  1/r  has  units  of  1/s. 


Polarization  Model  Likelihoods 

To  construct  the  likelihood  rates  p |_  and  p. ( for  the  polarization  model,  we  recall  that  the 

polarization  is  the  dipole  density  which  has  units  of  C/m2.  We  employ  the  energy  density  G  given 
by  (6)  which  corresponds  to  the  dipole  energy  (17).  Furthermore,  we  consider  a  continuum  of  dipole 
strengths  and  approximate  sums  by  integrals.  Finally,  we  approximate  evaluation  at  the  unstable 
equilibrium  Pm  by  evaluation  at  the  inflection  points  —Pi  or  Pj. 

These  approximations  yield  the  expressions 

fi(G)  =  Ce~GV/kT  (20) 


and 


P+-  = 


e~G(E,PI)V/kT 


7i 


l-G(E,P)V/kT  erfc x(Ep(t)) 


P-+  = 


l  e-G{E,-PI)V/kT 
"-Ple-G{E,P)V/kTdp 


7i 


erfc x(En(t)) 


(21) 


i  Pi 


for  the  Boltzmann  transition  probability3  and  transition  likelihood  rates.4  Here  erfcx  denotes  the 
complementary  error  function.  As  detailed  in  [13], 

En(t)  =  72 (E(t)  -  Ec )  ,  Ep(t)  =  72 (~E(t)  -  Ec), 

and 

_1  /  2Vp  _  1  [2  _  I  V  1  [ kT 

t  V  7 xkT  /3t\  ir  ’  y  2 kTr)  (5p\/2  ’  y  pV' 


Activation  Energies  A Ga 

Use  of  the  theory  of  thermally  activated  processes  to  construct  transition  likelihoods  has  the 
advantage  that  it  can  be  directly  motivated  by  fundamental  thermodynamics  concepts.  However, 
its  efficiency  diminishes  when  additionally  quantifying  90°  switching  or  2-D  or  3-D  polarizations 

3The  inclusion  of  the  reference  volume  V  in  the  expressions  (20)  and  (21)  is  due  do  the  definition  of  polarization  as 
the  dipole  density  and  G  as  an  energy  density  —  compare  with  the  dipole  relations  (18)  and  (19).  For  some  applications, 
the  estimation  of  parameters  yields  values  of  V  that  correspond  to  domain  or  grain  dimensions.  In  general,  however, 
it  is  simply  the  reference  volume  in  the  definition  of  P. 

4Because  probabilities  for  continuous  densities  are  defined  in  terms  of  integrals,  the  probability  of  measuring  a 
discrete  point  is  zero.  Hence  one  must  use  care  when  interpreting  the  relations  (20)  and  (21).  The  point  evaluations 
are  associated  with  the  discrete  set  of  polarization  states  with  sums  approximated  by  integrals  to  permit  formulation 
in  terms  of  error  functions.  The  proportionality  factor  associated  with  converting  sums  to  integrals  is  incorporated  in 
the  term  -.  An  alternative  interpretation  for  continuous  densities  is  provided  in  [75]. 
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due  to  the  higher  degree  of  energy  landscapes  and  difficulty  evaluating  integrals  based  on  inflection 
lines.  This  motivates  consideration  of  an  alternative  formulation  based  on  activation  energies  A Ga 
as  proposed  in  [47, 92] . 

As  illustrated  in  Figure  16,  the  activation  energy  is  the  difference  between  stable  and  unstable 
equilibria.  For  switching  between  the  negative  and  positive  well,  we  employ  the  formulation 


A  Ga_+{E) 


'  GU(E,  PM)  —  G-(E,  P~) 

<  Gu(-Ec,Pm)-G-(-Ec,P~) 
0 

f  \^(E-Ecf  ,  E<EC 
\  0  ,  E>  Ec 


—Ec  <  E  <  Ec 
E  <  -Ec 
E>  Ec 


where  GU,G_  are  defined  in  (6)  and  P~,Pm  are  defined  in  (7).  As  illustrated  in  Figure  16,  the 
unstable  equilibrium  exists  only  for  —  Ec  <  E  <  Ec.  For  computational  purposes,  we  employ  the 
limiting  values  for  E  >  Ec  and  E  <  —Ec.  This  differs  from  the  linear  formulation  employed  in  [47] 
but  yields  the  same  computational  value  for  the  transition  likelihoods. 

For  the  polarization  model  based  on  positive  and  negative  switching,  it  is  easy  to  specify  the 
unstable  component  GU(E,P)  which  provides  a  continuously  differentiable  transition  between  the 
stable  negative  and  positive  equilibrium  states.  For  2-D  and  3-D  polarization  models  or  models  that 
incorporate  90°  switching  behavior,  it  is  difficult  to  construct  continuously  differentiable,  or  even 
simply  continuous,  transitions  in  landscapes  between  the  regions  of  local  minima.  This  is  the  same 
problem  that  plagues  the  construction  of  multi-dimensional  splines.  Moreover,  it  is  unnecessary  since 
the  transition  regions  represent  unstable  dipole  behavior.  An  alternative  is  to  formulate  transition 
likelihoods  based  on  thermodynamic  driving  forces. 

For  this  regime,  the  thermodynamic  driving  force  is  defined  to  be 


F_+(E)  =  G-(E)  -  G+(E)  =  2EPr 
F~\ —  (E)  =  G+(E)  -  G-(E)  =  -2EPr 


(see  [9]  and  included  references).  If  we  define  the  critical  driving  force  to  be 


Fc  =  2  EcPr, 


(22) 


then  the  activation  energy  can  be  expressed  as 


A  Ga_+{E) 


AGo  (1  —  F_+(E)/FC)2 
0 


F_+(E)  <  Fc 
F>FC 


where  A  Go  =  Fc/4  denotes  the  value  of  the  energy  barrier  at  zero  driving  force.  The  definition  of 
AG“  _  is  analogous.  We  note  that  in  the  homogenized  energy  model,  Fc,  and  hence  AGo,  is  treated 
as  a  material  parameter  whose  values  are  realizations  of  an  underlying  density. 

The  likelihood  relations 

p_+(E)  =  -e~AG-+^v/kT  ,  p+(E)  =  -e~AGa+-^v/kT  (23) 

r  r 

follow  directly  from  the  previous  Boltzmann  theory  when  we  consider  only  polarization  values  P~,  P+ 
and  Pm-  This  formulation  facilitates  implementation  by  eliminating  integration  over  complex  regions 
but  it  neglects  the  geometry  of  the  energy  surface  in  regions  near  minima  since  it  is  equivalent  to 
summing  or  integrating  over  square  energy  wells.  In  addition  to  being  used  in  [47],  the  activation 
energy  is  employed  in  [4]  to  construct  Debye  relations  for  ferroelectric  materials.  It  is  also  used  in  [10] 
for  constructing  magnetic  models  and  [20]  where  analogous  models  are  developed  for  SMA. 
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3.3  Polarization  Kernel  that  Incorporates  Thermal  Relaxation 

Solution  of  the  differential  equation  (10)  using  the  likelihood  relations  (21)  or  (23)  quantifies  the 
dipole  fraction  x+  and  x_  =  1  —  x+.  The  general  polarization  hysteron  is  then  given  by  (11). 

As  detailed  in  [75,  78-80],  P~  and  P+  for  the  thermally  active  case  are  given  by 


P+(E)=  /  Pp(G(E,  P)dP 
Jpi 


— P , 


P~(E)  = 


P/i(G(E,  P)dP 


(24) 


where  //  is  defined  in  (20).  For  the  likelihood  relations  (21),  it  is  shown  in  [13]  that  the  kernel  can 
be  expressed  as 

P(E;  Ec)=E-Pr  +  2  x+PR  +  P(E) 

V 

where  P(E ;  Ec)  =  74  (x+  —  1)  p _ +  p-\ _ x+  and  74  =  To  determine  the  limiting  behavior  of  P, 

we  note  that  when  p _ |_  is  large  and  hence  p. \ _ is  small,  the  likelihood  that  dipoles  are  positive  is  large 

which  implies  that  x+  ~  1.  Similarly,  p [_  small  and  p. \ large  yields  x+  ~  0.  Finally,  P(0)  =  —74 

where  74  is  typically  on  the  order  of  ICR6  to  ICR8  for  ferroelectric  materials.  The  approximation 
P(E)  =  0  then  yields  the  kernel  expression 


P(E,  Ec )  = - Pr  +  2 x+Pr 

V 


(25) 


which  is  the  negligible  thermal  activation  relation  (13).  Alternatively,  it  is  shown  in  Section  2.6.3 
of  [75]  that  use  of  Dirac  sequences  yields  (25)  as  an  appropriate  low  thermal  activation  limit. 


3.4  Homogenized  Energy  Model 

To  incorporate  the  effects  of  polycrystallinity,  material  nonhomogeneities,  and  variable  interaction 
fields,  we  assume  that  interaction  fields  Ej  and  certain  material  coefficients  are  manifestations  of 
underlying  densities  rather  than  constants.  It  is  illustrated  in  [13,75,78,80]  that  the  assumption  that 
coercive  fields  are  distributed  yields  a  ±180°  switching  polarization  model 


p(E(ty,x°+)  = 


Jo 

E 


P(E(t)  +  Ep,  Ec)vI(EI)i'c(Ec)dEIdEc 


(26) 


—  Pr  +  2  Pr  / 

V  Jo 


x+  (Ee (t ) ;  Ec)  V! (Ei )  vc (Ec) dEj dEc , 


that  accurately  quantifies  polycrystalline  behavior.  Here 

Ee(t)  =  E(t )  +  Ej 


(27) 


is  the  effective  electric  field,  is  the  initial  fraction  of  positively  oriented  dipoles,  and  vc  and  uj  are 
densities  associated  with  the  coercive  and  interaction  fields  which  satisfy  the  constraints 


(i)  vc(Ec)  defined  for  Ec  >  0, 

(ii)  vi(-Ei)  =  pj(Ei), 

(hi)  Wc(Ec)\  <  c1e~aiEc  ,  \pi{EI)\<c2e~a^El^ 


(28) 


for  positive  constants  ci,  04,02,02-  These  assumptions  enforce  the  physical  properties  that  local 
coercive  fields  are  positive,  low-field  Rayleigh  loops  are  symmetric,  and  local  coercive  and  interac¬ 
tion  fields  decay  as  a  function  of  distance.  Details  regarding  density  construction  are  provided  in 
Section  4.2.2  and  quadrature  rules  to  approximate  the  integrals  in  (26)  are  provided  in  Section  4.3. 
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Figure  17:  Hysteresis  behavior  when  (a)  vc  is  distributed  but  vj  is  constant,  and  (b)  vj  is  distributed 
but  vc  is  constant. 


The  necessity  of  considering  densities  for  both  Ej  and  Ec  is  illustrated  in  Figure  17  by  illustrating 
the  hysteresis  behavior  that  results  if  either  is  neglected.  As  shown  in  Figure  17(a),  the  assumption 
that  vc  is  distributed  but  vj  is  constant  yields  loops  which  exhibit  no  switching  as  fields  are  reduced 
until  the  value  E  =  0  with  similar  behavior  for  increasing  fields.  As  illustrated  in  Figure  ll(d)-(f) 
for  PZT,  this  is  not  the  case  for  ferroelectric  materials.  The  assumption  that  vc  is  constant  and  vi 
distributed  yields  the  minor  loop  behavior  depicted  in  Figure  17(b).  Due  to  the  symmetry  of  vj,  this 
forces  minor  loops  to  be  rotationally  symmetric  about  E  =  0  which  is  not  the  behavior  exhibited  by 
materials.  This  necessitates  the  inclusion  of  both  densities. 


4  Homogenized  Energy  Model  —  Strain  and  Polarization 

To  quantify  the  strain  and  polarization  behavior  detailed  in  Section  2  for  various  stress  regimes,  it  is 
necessary  to  additionally  incorporate  90°  switching.  As  motivated  by  common  actuator  geometries 
and  to  simplify  the  discussion,  we  still  focus  on  1-D  input  and  output  responses.  For  tetragonal 
materials,  the  projection  of  P  in  Figure  6  onto  the  3-axis  yields  spontaneous  polarization  variants 
or  states  at  ±180°,  90°  which  are  respectively  designated  P^,a  =  ±,90.  Similarly  rhombohedral 
materials,  such  as  PLZT,  have  variants  at  ±180°,  70.5°  and  109.5°. 

We  note  that  the  construction  of  grain-level  or  single  crystal  polarization  and  strain  kernels 
follows  the  development  in  [47,92],  However,  these  authors  do  not  consider  the  homogenized  energy 
framework,  developed  in  Section  4.2,  which  is  necessary  to  characterize  polycrystalline  materials. 


4.1  Strain  and  Polarization  Kernels 


For  a  =  ±180,  90,  the  Helmholtz  and  Gibbs  energy  densities  for  the  a- well  are 


i>a(P,£) 


2 vi(P  -  Pr ?  +  -Ya(c  -  Cr)2  +  K(P  -  P%)(e  -  e°R) 


(29) 


and 

Ga(E,a;  P,e)  =  ipa(P,£)  —  EP  —  ae.  (30) 

As  summarized  in  the  nomenclature  table  at  the  beginning  of  the  paper,  PR  and  are  the  remanent 
polarization  and  strain  of  the  a- variant.  The  parameters  77^,  and  ha  are  the  inverse  susceptibility 
at  constant  strain,  elastic  stiffness  at  constant  polarization,  and  piezoelectric  constant,  respectively. 
For  a  fixed  applied  stress  a  and  field  E,  the  conditions 


dG  dG 
HP  _  1  ’  ~de  ~ 
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can  be  used  to  determine  pairs  (P,n,£m)  that  minimize  the  Gibbs  energy.  Specifically,  this  yields 


E  =  rfa(P%-Pg)  +  ha(e%-e%) 
o  =  Y?{sam  -eaR)  +  ha(P«-P%) 


which  can  be  inverted  to  obtain 


Pm  —  Pr  +  XqP  +  da° 
£m  =  £r  +  daE  +  sp  a 


(31) 


where 


YP 


Xa 


,  da  — 


hn 


To 


’  yaPV£a  -  K 


Ya  'nEa  -  hi  ’  III  -  Y£rfa 
Note  that  (31)  provides  the  linear  domain-level  constitutive  equations  shown  in  Figure  4.  The 
minimum  of  the  Gibbs  energy  density  in  each  a-well  can  thus  be  expressed  as 


Gam(E,  a)  =  -\xaaE2  -  \sp<?2  -  daEa  -  EP%  -  aeaR. 


(32) 


We  note  that  within  each  well,  the  Gibbs  energy  density  (32)  is  the  negative  Legendre  transform  of 
the  Helmholtz  energy  density  (29). 

To  construct  the  activation  energies  employed  in  transition  likelihoods,  we  note  that  the  thermo¬ 
dynamic  driving  forces  required  to  transition  from  an  a-well  to  /3-well  can  be  expressed  as 


Fap{E,a)  =  Gam(E,  a)  —  Ggm(E,  a ) 

=  lAx'u-E2  -  WfX  -  A dagEa  -  EApf  -  aAef 

where 

AXa/3  =  xl  -  xh  ,  A  sf„  =  sf  -  sf  .  A  da0  =  d<,-dn 

A pf  =  p%-pL  A4”  =  4  -  4- 

Based  on  the  assumption  that  dipoles  are  restricted  to  three  orientations  a  =  ±,  90,  we  let  x+,  X- 
and  xgo  denote  the  dipole  fractions  associated  with  negatively,  positively,  and  90°-oriented  dipoles. 
The  evolution  of  dipole  fractions  is  governed  by  the  differential  equation 


x-  =  ~(p~  90  +  P-+)x _  +  pgo-.xgo  +  P+-X+ 
X90  =  P-90X-  -  (P9 0-  +  P90+)®90  +  P+90%+ 
x+  =  P-+X-  +  P90+X90  -  {p+ 90  +  p+-)x+ 


which  can  be  simplified  to 


X-  =  -(P-90  +P~+  +P90-)X-  +  (p+-  ~P90-)x+  +P90- 
x+  =  (p-+  ~P90+)X-  -  (P+90+P+-  +  P90+)x+  +  P90+ 


(33) 


using  the  identity  x+  +  X-  +  X90  =  1. 

As  motivated  by  the  discussion  in  Section  3.1,  the  likelihood  rates  pap  of  transitioning  from  an 
a-well  to  /3-well  are  specified  by 


Pap{E,c r) 


1  c-AG*„(E,o-)V/kT ' 
Ta(3 


(34) 
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The  activation  energy  is  specified  by  the  relation 


A  Gaap(E,a;Fc) 


f  AG0(1  —  Fag(E,a)/Fc)2 

l  0 


Fa/3(E,a)  <  Fc 
Fag(E,a )  >  Fc 


,  AG0  = 


\Fc  ,  180°  Switching 
YqFc  ,  90°  Switching 


where  A  Go  is  the  energy  barrier  at  zero  driving  force.  We  note  that  the  critical  driving  force  Fc  is 
assumed  to  be  a  manifestation  of  an  underlying  density  in  the  homogenized  energy  model. 

It  is  noted  in  Section  2.2  and  [47,96]  that  whereas  180°  switchings  is  often  comprised  of  successive 
90°  switches,  it  typically  occurs  on  a  faster  timescale  than  strain-producing  90°  switches.  These 

different  timescales  are  accommodated  by  the  inclusion  of  180°  likelihood  rates  p _ _ and  use  of 

different  relaxation  times  Ta$  for  90°  and  180°  switchings. 

The  polarization  and  strain  kernels  are  given  by 


p=  Y  x<*pa 

a=±,9  0 

where 

Pa  =  f  P/i(P,  e)dPde 

J  a 

are  the  polarization  and  strain  associated  with  each  well.  In  (35),  integration  is  performed  over 
a  region  of  the  energy  minima  and  p(Pa,£a)  given  by  (20)  quantifies  the  probability  of  finding  a 
specific  polarization  and  strain  pair  (Pa,£a). 

The  evaluation  of  the  integrals  (35)  is  difficult  for  higher-dimensional  energy  landscapes  which 
significantly  diminishes  implementation  speeds.  To  address  this,  we  employ  the  low  thermal  acti¬ 
vation  approximations  discussed  in  Section  3.2  and  employ  the  approximations  Pa  =  P^,£a  = 
where  are  given  in  (31).  The  kernels  are  subsequently  given  by5 


,  £  =  Y  Xa£C 

a=±,9  0 


,  £a  =  /  ep(P,  e)dPde 


(35) 


p=  Y  x«p™ 

a=±,90 


£  = 


£ 

a=±,90 


Xn£r 


(36) 


To  construct  a  form  of  the  kernel  that  facilitates  comparison  with  existing  models  posed  in  terms 
of  reversible  and  irreversible  components  and  construction  of  constitutive  models  for  coupled  and 
distributed  structures,  we  consider  (36)  in  light  of  the  assumption  that 


X+  =  X-  =  X90  =  X 

S+  =  S-  =  S90  =  sE 

Pr  =  0)  Pr  =  ~Pr  j  eR  =  eR, 

dgo  =  0,  d-  =  -d+ 


(37) 


790-  —  r-90  —  790+  —  7+90  —  7go  ,  T__+  —  T+_  —  TigO- 


This  assumption  is  based  on  the  observed  material  behavior  and  does  not  reduce  the  model’s  gener¬ 
ality.  Substitution  of  (37)  into  (36)  yields  the  relations 


P(E,  a)  =  d(E ,  a)a  +  xaE  +  Pirr(E,  a) 
e(E,  a)  =  sEcr  +  d(E ,  a)E  +  £irr(E,  a) 

5  The  input  and  model  parameters  exhibit  the  dependencies 

E(t),  a(t),  x±,90(t;  E,  a ,  Fc),pap{E ,  a;  Fc),  P(E,  a;  Fc),s(E,  a\  Fc). 
Where  the  meaning  is  clear,  we  suppress  certain  dependencies  to  simplify  subsequent  discussion. 


(38) 
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where 


Pirr(E,a)  =  ^  P%xa{E,a ) 

a=±,90 

£irr{E,(j)  =  'y  '  £j^Xa{E1  O')  (39) 

o=±,90 

d(E,a)=  ^2  daxa(E,cr) 
a=±,9  0 

are  the  grain-level  polarization,  strain  and  coupling  relations  illustrated  in  Figure  4. 

Example  1. 

The  kernel  relations  (39)  can  be  used  to  model  the  behavior  of  certain  single  crystal  compounds 
if  interaction  fields  are  negligible.  This  is  illustrated  in  Figure  18  to  demonstrate  the  ferroelastic 
switching  behavior  discussed  in  Section  2  and  depicted  in  Figure  9.  The  parameters  in  Table  1  are 


(a) 


(d) 


(g) 


_0.6 

^0.4 

2 

55  0.2 


0 


-1.5  -1  -0.5  0  0.5  1  1.5 

E  (MV/m) 

(b) 


(e) 


(h) 


E  (MV/m) 


(c) 


(f) 


-0.4  0  0.4 

P  (C2/m) 


(i) 


0)  (k)  (£) 

Figure  18:  Simulated  polarization  and  strain  for  fixed  stress  levels  of  (a)-(c)  0  MPa,  (d)-(f)  -4  MPa, 
(g)-(i)  -8  MPa,  and  (j)-(-Z)  -30  MPa  with  instantaneous  switching  of  xa. 
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p+ 

R 

4 

P90 

eR 

X+ 

d+ 

4 

t90 

T1S0 

0.26 

6.7  x  10"3 

-4.2  x  10“3 

5.0  x  10"8 

374  x  10“12 

18.8  x  10"12 

1.0  x  10“3 

3.0  x  10-3 

Table  1:  Parameters  used  in  the  single  crystal  simulations. 


(m)  (n)  (o)  (p) 

Figure  19:  Gibbs  energy  density  for  fixed  stress  levels  of  (a)-(d)  0  MPa,  (e)-(h)  -4  MPa,  (i)-(£) 
-8  MPa,  and  (m)-(p)  -30  MPa  at  the  indicated  input  field  levels. 


representative  of  those  for  standard  ferroelectric  materials.  The  corresponding  Gibbs  energy  density 
landscapes  are  plotted  in  Figure  19.  It  is  observed  that  the  application  of  increasing  compressive 
prestresses  shifts  the  landscape  to  favor  90°  initial  orientations. 


4.2  Homogenized  Energy  Model 


As  noted  in  Section  3.3,  the  effects  of  polycrystallinity  and  variable  interaction  fields  can  be  incor¬ 
porated  by  considering  certain  parameters  to  be  manifestations  of  underlying  densities.  Here  we 
employ  the  homogenized  energy  formulations 


P(E(t),a(t);x°+) 


roc  roc 

/  /  -p(Ee(t),cr(t);  Fc)vI(EI)i,c(Fc)dEIdFc 

J  0  J  —  oo 


'  0  J  —  oo 

roc  roc 

1 0  J — oo 


e{Ee(t),(r(t)\Fc)vI(EI)vc(Fc)dEIdFc 


(40) 
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where  the  effective  field  Ee  is  defined  in  (27).  The  choice  of  the  critical  driving  force  density  vc(Fc ) 
is  motivated  by  its  role  when  constructing  activation  energies  for  likelihood  construction  and  its 
relation  (22)  to  the  coercive  field  Ec.  The  densities  satisfy  the  criteria  (i)-(iii)  in  (28)  and  their 
construction  is  discussed  in  Section  4.2.2. 

With  the  assumption  (37)  and  resulting  kernel  relations  (38),  (40)  can  be  expressed  as 

P{E,  a)  =  d(E,  <r)a  +  x^E  +  Pirr{E,  u) 

_  (41) 

e(E,  a)  =  sh a  +  d(E,  a)E  +  £irr(E,  cr) 


where 

poo  poo 

d{E,  a )  =  /  /  d(Ee]Fc)uI(EI)iyc(Fc)dEIdFc 

J  0  J — oo 

poo  poo 

Pirr{E1(j)=  /  /  Pirr(Ee-,Fc)uI(EI)uc(Fc)dEIdFc 

J  0  J  —oo 

poo  poo 

£irr{E^(J^  —  /  I  £irr(Ee]  Fc')l/i(Ej')l/c(FcS)dEidFc. 

J  0  J  —oo 

We  note  that  the  terms  xaE  and  sEa  result  from  the  fact  that 


uI(EI)vc(Fc)dEIdFc  =  1. 


(42) 


Techniques  to  evaluate  the  integrals  in  (42)  are  detailed  in  Section  4.3. 

We  note  that  (41)  has  the  form  of  the  general  constitutive  relations  (3)  that  others  have  derived 
using  various  micromechanical  and  phenomenological  theories;  see  also  Figure  4.  The  homoge¬ 
nized  energy  model  provides  the  advantage  of  nricronrechanical  energy  analysis  in  combination  with 
stochastic  homogenization  techniques  to  facilitate  efficient  model  implementation  and  data-driven 
parameter  estimation  techniques  [31]. 


4.2.1  Determination  of  Initial  Hysteresis  States 

Implementation  of  the  model  (40)  requires  the  specification  of  the  initial  dipole  states.  Due  to  the 
multivalued  nature  of  hysteresis  phenomena,  this  a  shared  requirement  for  all  hysteresis  models. 
This  necessitates  determining  the  initial  state  of  devices  with  hysteretic  components  so  it  is  an 
experimental  and  implementation  issue  in  addition  to  being  a  modeling  requirement. 

One  technique  to  address  this  is  to  experimentally  drive  devices  to  positive  or  negative  saturation 
before  the  specified  task  to  ensure  an  initial  remanence  state.  To  initialize  the  model,  one  can  start 
with  arbitrary  initial  dipole  fractions;  e.g.,  x^_  =  Xg0  =  x^_  =  |.  Driving  to  positive  saturation  yields 
x+  =  1 ,  Xg0  =  x^_  =  0  and  returning  the  field  to  zero  yields  remanence  dipole  fractions  xr+ ,  Xg0  and 
xr_  that  can  be  used  as  initial  values  for  subsequent  characterization  and  simulations. 


4.2.2  Density  Construction 

A  critical  issue  when  constructing  models  for  specific  applications  concerns  the  construction  of  the 
densities  vi(Ei)  and  vc(Fc ).6  Here  we  discuss  a  method  introduced  in  [24,30]  that  employs  the 

6As  illustrated  in  Figure  11,  the  polarization  burst  region  and  strain  elbow  occur  at  basically  the  same  field  Ec. 
This  indicates  that  one  value  of  Fc  and  density  ^(A)  can  be  used  for  both  180°  and  90°  switching  processes. 
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expansions 


Ka 


Vc(Fc)  —  ^  '  Q7,0fc(-Fc) 

r  i  z — ' 


Cl 

1 


k=ka 

Ka 


(43) 


vi(Ei)  =  -  V  (3kVk{Ei) 
c2  ' 


k=kp 


where  the  basis  functions  4>k(Fc )  and  ipk(Ej)  are  lognormal  and  normal  functions  and  the  coefficients 
CKfc  and  (3k  are  determined  through  a  least  squares  fit  to  data  as  detailed  in  [31].'  The  constants 


A' 


Ka 


Cl  —  ^2  ak  1  C 2  —  ^  /3fc 


k — fen 


k=ka 


ensure  integration  to  unity. 

The  interaction  field  basis  functions  are  defined  by 


Vk(Ej)  = 


1 


Gk1\f7lTX 


-E(/2{akI)'2 


where  the  standard  deviations  are  taken  to  be 

07  =  {o/fc}  =  {2*07}  j  k  =  ka,  ■  ■  ■ ,  Ka. 

As  detailed  in  [30,31],  data-driven  techniques  can  be  used  to  obtain  initial  values  of  07. 
Since  Fc>  0,  we  employ  lognormal  functions 


(44) 


4>k{Fc )  = 


1 


B—  HR)- Mc]2/2Kfc)2 


akFcy/2TT 


for  the  driving  force  density.  Data-driven  techniques  to  determine  /rc  =  lnFc  are  detailed  in  [31]. 
The  standard  deviations  ac  are  specified  in  a  manner  analogous  to  (44).  We  note  that  for  typical 
applications,  reasonable  accuracy  can  be  obtained  with  ka  =  kg  =  —3,  Ka  =  Kg  =  1.  Representative 
basis  functions  are  plotted  in  Figure  20. 

7The  expansions  (43)  are  linear  with  respect  to  the  parameters  Ofe  and  /3k-  This  permits  the  use  of  linear  least  squares 
and  linear  adaptive  techniques  for  parameter  estimation  and  adaptation.  However,  these  techniques  are  constrained 
by  the  requirement  that  parameters  are  nonnegative. 


Figure  20:  (a)  Normal  basis  functions  (pk(Ej)  and  (b)  lognormal  functions  (f>k(Fc)  used  in  the  expan¬ 
sions  (43)  for  vj(Ej)  and  uc(Fc). 


32 


4.3  Quadrature  Techniques 

To  implement  the  model,  the  integrals  defining  the  terms  d,  Pirr  and  £irr  in  (42)  must  be  approxi¬ 
mated  in  a  manner  that  is  both  efficient  and  accurate.  We  first  note  that  by  employing  the  relations 
(39),  these  terms  can  be  expressed  as 


Pirr(E,  <7 ) 

Sirr(E  ■  a) 
d(E,  a)  = 


=  E 

a=±,90 


xa[Ee\ Fc)uI(Ei)i/c(Fc)dEjdFc 


-  E 

a=±,90 


xa{Ee ;  Fc)uI(EI)fc(Fc)dEIdFc 


^  dQ  /  xa(Ee;Fc)uI(EI)iyc(Fc)dEIdFc 

q=±,90  •'° 


so  the  integration  involves  the  dipole  fractions  in  a  manner  analogous  to  (26).  Secondly,  the  densities 
satisfy  the  exponential  decay  constraints  (28)  as  shown  in  Figure  20.  Hence  approximation  algorithms 
can  be  defined  on  finite  domains  rather  than  necessitating  quadrature  techniques  for  infinite  and 
semi-infinite  domains.  We  illustrate  here  a  midpoint  formula  and  refer  the  reader  to  [75]  for  details 
regarding  composite  Gaussian  quadrature  techniques.  As  detailed  in  [31],  truncation  of  the  domains 
and  use  of  composite  techniques  proves  more  efficient  to  implement  than  Gauss-Hermite  and  Gauss- 
Laguerre  algorithms,  that  are  designed  specifically  for  infinite  and  semi-infinite  intervals,  because 
they  allow  accurate  evaluation  of  the  density  expansions  (43)  and  construction  of  lookup  tables  to 
improve  algorithm  efficiency. 

If  we  let  FCi,Ejj  and  Vi,Wj  respectively  denote  the  quadrature  points  and  weights,  we  obtain  the 
approximate  relations 


Ni  N, 

xa(E  +  Ej]  Fc)zv/(H7)z/c(Fc)dT/dFc  «  EE-  a{E  +  ETj  ;  FCi)uI( Ej.  )vc(FCj)viWj . 

i=  1  j= 1 

For  the  validation  results  reported  in  Section  6,  sufficient  accuracy  was  achieved  using  the  trapezoid 
rule  with  Nt  =  Nj  =  41. 


5  Constitutive  Relations  for  Structural  Models 

For  isolated  materials,  the  polarization  and  strain  relations  (40)  or  (41)  may  be  suitable  for  character¬ 
ization  and  simulations.  However,  for  most  applications,  the  ferroelectric  materials  are  components 
in  a  system  which  necessitates  the  development  of  system  models  employing  (40)  or  (41).  Here 
we  discuss  the  construction  of  a  simple  lumped  model  for  PZT  employed  with  a  spring,  damping 
element,  and  prestress  and  the  development  of  constitutive  equations  that  can  be  used  to  construct 
rod,  beam,  shell  or  membrane  models  for  systems  with  ferroelectric  actuators. 

5.1  Lumped  Actuator  Model 

We  first  consider  a  PZT  actuator  subjected  to  a  prestress  <to  and  restoring  spring  with  stiffness  k 
as  illustrated  in  Figure  21(a).  This  is  representative  of  characterization  experiments  reported  in  [92] 
and  can  provide  a  simplified  model  for  certain  applications.  The  actuator  is  assumed  to  have  length 
L.  cross-sectional  area  A,  and  longitudinal  displacements  of  the  actuator  end  are  denoted  by  u(t). 
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For  this  configuration,  the  inputs  or  loads  are  the  field  E  and  prestress  tq.  Hence  the  Gibbs 
energy  density  (30)  or  (32)  is  the  appropriate  thermodynamic  potential.  The  results  strains  are 
specified  by  (40)  or  (41). 

The  total  applied  force  is 

F  =  FS  +  F0  (45) 


where 


Fs  =  —ku  ,  Fq  =  o-qA  ,  with  a q  <  0, 


are  the  forces  due  to  the  restoring  spring  and  prestress.  Since  e  =  it  follows  that  the  applied 
stress  is 

(T=  —^-k  +  a0.  (46) 

The  substitution  of  (46)  into  (41)  yields  the  relation 


A 

A  +  sELk 


sE  tq  T  d  (-F,  — Lek/ A  up)  E  -f-  £  (E ,  — Lsk  j A  uq) 


(47) 


quantifying  strains  in  terms  of  E  and  cto-  Due  to  its  implicit  nature,  solution  of  (47)  requires  a 
nonlinear  iterative  method  such  as  a  fixed-point  or  Newton-based  algorithm. 

For  many  applications,  however,  the  spring  stiffness  k  is  sufficiently  small  that  |  k |  <<  | cxq | - 
For  these  cases,  (47)  can  be  approximated  by  the  explicit  relation 

£  =  4  +AgE Lk  [^O  +  d  (E'  -0)  E  +  £irr(E ,  CT0)]  .  (48) 


It  is  illustrated  in  Section  6.1  that  for  data  from  York,  spring  constants  on  the  order  of  10  times 
physical  values  are  required  before  (48)  begins  to  differ  significantly  from  (47).  For  regimes  that 
require  (47),  the  iterative  fixed-point  algorithm 


£n+1  =  4  +AsELk  [sE°o  +  d  (E, -L£nk/A  +  £j0)  E  +  £  (E,  -L£nk/A  +  a0)] 

can  be  used  to  efficiently  compute  strains  as  a  function  of  E  and  <ro- 

To  incorporate  the  damping  element  depicted  in  Figure  21(b),  we  include  a  damping  force  F(j  = 
—cu  in  (45)  to  obtain  the  applied  stress  relation 


— Le  Li 

<J  =  — T—k - —c  ■ 

A  A 


cr0- 


This  then  yields  the  differential  equation 

seLc  .  A 

-£  H-  £  — 


A  +  sELk  A  +  sELk 

which  specifies  strains  produced  by  the  actuator. 


sEa0  +  d(E,ao)  E  +  £irr  {E,<J  o)] 


(49) 


(a)  (b) 


Figure  21:  (a)  PZT  actuator  with  a  prestress  uo  and  restoring  spring,  and  (b)  actuator  that  addi¬ 
tionally  incorporates  a  damping  element. 
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5.2  Constitutive  Relations  for  Distributed  Models 

The  nature  of  inputs  differs  somewhat  for  distributed  structures  such  as  rods,  beams,  shells  and 
membranes.  To  illustrate,  consider  a  rod  of  length  L  as  illustrated  in  Figure  22.  The  applied  loads 
again  are  the  field  E  and  prestress  <7o  at  the  end  of  the  rod.  The  development  of  PDE  structural 
models  requires  the  quantification  of  stresses  along  the  rod  length. 

To  determine  whether  strains  or  stresses  should  be  considered  as  inputs  for  0  <  x  <  L,  we 
consider  connections  at  x  =  a  that  block  rod  movement.  The  stress  do  produces  strains  e(t,x)  and 
displacements  u(t,x )  for  x  >  a  along  with  a  stress  cro(t,a).  Despite  this  stress,  there  are  no  strains 
nor  displacements  for  x  <  a  due  to  the  blocking  elements.  Hence  strains  rather  than  stresses  should 
be  considered  as  inputs  along  the  rod  so  the  independent  variables  for  0  <  x  <  L  are  E  and  e. 

For  this  independent  variable  set,  the  appropriate  thermodynamic  potential  is  the  electric  Gibbs 
energy  density 

Gea  =  ifia  —  EP  =  Ga  +  ere 

(see  pages  63-64  of  [75]).  To  obtain  appropriate  constitutive  relations,  one  can  repeat  the  analysis 
of  Section  4.1  using  G®  rather  than  Ga.  While  fundamental  in  nature,  this  approach  obscures  a 
comparison  with  the  constitutive  relations  (41)  used  to  model  (E,ao)  inputs.  Instead,  we  directly 
reformulate  the  relation  (41)  to  provide  stress  as  a  function  of  strain. 

To  obtain  an  explicit  constitutive  equation,  we  assume  that  |<7o|  >>  |cr|  for  d  and  £irr  defined  in 
(42)  and  employ  d(E,ao)  and  £irr(E,aQ).  Inversion  of  the  strain  relation  in  (41)  then  yields 

cr(E,  £)  =  YE£  -  e(E,  a0)E  -  YE£irr(E ,  a0)  (50) 

where 

YE  =  Je  ’  e(E,a0)  =  ^d(E,a0). 

The  constitutive  relation  (50)  can  be  used  to  construct  rod,  beam  and  shell  models  for  distributed 
structures.  To  illustrate,  consider  a  rod  having  density  p  and  cross-sectional  area  A.  Longitudinal 
displacements  are  denoted  by  u(t,x).  As  detailed  in  Chapter  7  of  [75],  force  balancing  yields 

d2u  dN 
P  dt 2  dx 

where  the  force  resultant  is  N  =  a  A.  The  stress  is  given  by  (50). 

For  general  electro-elastic  structures,  one  would  employ  the  relations 

V  •  D  =  0  ,  D  =  eo  E  +  P 

pu  =  V  •  cr 

VxB  =  0  ,  E  =  —  V</?. 

The  use  of  (50)  to  construct  models  for  PZT-based  macro-fiber  composite  (MFC)  actuators  operating 
in  highly  hysteretic  and  nonlinear  regimes  is  detailed  in  [32], 


L _ 

J _ 

x^a  u(t,x) 
(a) 


L _ _ 

_J _ -::O0 

x=a  u(t,x) 

(b) 


Figure  22:  (a)  Rod  of  length  L  with  a  connection  at  x  =  a  that  blocks  movement,  and  (b)  displace¬ 
ment  produced  by  a  prestress  cq- 
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6  Model  Validation 


The  relations  (40)  or  (41)  characterize  the  polarizations  and  strains  generated  by  fields  or  stresses 
applied  to  isolated  ferroelectric  materials.  For  actuators  subjected  to  a  prestress  and  restoring  spring, 
the  relations  (47),  (48)  or  (49)  quantify  the  generated  stresses  whereas  (50)  provides  a  nonlinear  con¬ 
stitutive  relation  for  distributed  structures.  In  this  section,  we  illustrate  the  experimental  validation 
of  the  model  for  a  prestressed  PZT  element  with  a  restoring  spring,  rate-dependent  PZT,  and  PLZT 
subjected  to  various  compressive  prestresses. 

Details  regarding  highly  efficient  implementation  algorithms  and  data-driven  parameter  estima¬ 
tion  techniques  are  provided  in  [31] .  Additional  examples  illustrating  the  experimental  validation  of 
the  model  for  PZT  with  variable  loading  rates  and  single  crystal  BaTiC>3  are  provided  in  that  paper. 

6.1  Model  Validation  for  a  Prestressed  PZT  Actuator 

We  illustrate  here  validation  of  the  model  using  data  reported  in  [92]  for  a  prestressed  PZT  actuator 
with  a  restoring  spring.  The  multilayer  actuator  (NEC  Model  AE0505D08)  employed  the  soft 
PZT  material  PIC151  and  was  comprised  of  80  active  layers,  each  having  thickness  0.1  mm.  The 
dimensions,  including  packaging,  were  6.5  x  6.5  x  10  mm  which  yielded  the  cross-sectional  area 
A  =  42.25  mm2.  Displacements  were  measured  using  a  fiber-optic  displacement  sensor  having  an 
effective  resolution  of  0.1  /rm  and  strains  were  computed  by  dividing  the  total  thickness  (8  mm)  of  the 
active  layers.  The  data  used  for  model  validation  was  collected  with  a  prestress  of  <7o  =  —10.6  MPa 
with  a  restoring  spring  having  a  stiffness  k  =  2.7  N//mi.  The  polarization  was  measured  using  a 
Sawyer- Tower  circuit  comprised  of  a  reference  capacitor  connected  in  series  with  the  PZT  actuator. 
Strain  and  polarization  data  collected  at  a  loading  rate  of  0.5  AV-  is  plotted  in  Figure  23. 
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Figure  23:  Fit  of  the  polarization  and  strain  model  to  data  from  [92]  with  a  prestress  of  ao  = 

—10.6  MPa:  (a)  time  domain  and  (b)-(d)  phase  space;  ( - )  experimental  data,  (—•  —  •)  model 

fit,  ( - )  model  fits  for  8  selected  loading  cycles. 
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Figure  24:  (a)  Critical  driving  force  density  vc{F^) 
equally  spaced  quadrature  points  marked  as  dots. 
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and  (b)  interaction  field  density  vi(Ei)  with  41 


In  Figure  23(a),  it  is  observed  that  following  the  application  of  a  saturating  1.5  kV/mm  field,  fields 
were  decreased  and  held  at  15  different  levels  before  being  increased  back  to  saturation.  Depending 
on  the  held  field  value,  this  produced  varying  degrees  of  90°  switching,  that  is  manifested  in  the 
strain  measurements,  along  with  slower  creep  behavior  in  both  the  polarization  and  strain. 

The  polarization  was  modeled  by  the  relation  (41)  and  the  strain  by  (48)  with  the  measured 
field  plotted  in  Figure  23(a)  used  as  input.  We  neglected  damping  forces  and  the  use  of  (49)  due  to 
the  slow  loading  rate.  For  maximum  strains  on  the  order  of  emax  =  0.0014,  j  |  =  0.89  MPa 

<<  1 |  =  10.6  MPa,  thus  motivating  the  use  of  the  linear  relation  (48)  rather  than  the  nonlinear 
relation  (47).  Further  discussion  detailing  the  validity  of  the  linear  strain  relation  for  this  actuator 
model  is  provided  at  the  end  of  this  section. 

Use  of  the  data-driven  parameter  estimation  techniques  discussed  in  [31]  yielded  the  initial  and 
optimized  parameters  summarized  in  Table  2  and  the  resulting  model  fits  plotted  in  Figure  23.  From 
Figure  23(a),  it  is  observed  that  the  model  accurately  characterizes  both  the  switching  and  creep 
behavior  of  the  polarization  and  strain  as  a  function  of  time.  To  illustrate  the  E-P,  E-e  and  P-e 
behavior,  the  loops  indicated  in  Figure  23(a)  were  selected  and  plotted  in  Figures  23(b)-(d).  The 
estimated  interaction  and  coercive  field  densities  vi  and  vc  are  plotted  in  Figure  24.  The  reader  is 
referred  to  [31]  for  statistical  properties  of  the  densities. 
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Table  2:  Initial  and  optimized  parameter  values  for  the  90°  model  with  do  =  —10.6  MPa.  The 
optimized  density  parameter  values  aq  =  2.58,  «o  =  1.66,  a_i  =  2.01,  a_2  =  2.10,  a_3  =  2.03  and 
Pi  =  1.12,  /?o  =  0.71,  /?_i  =  2.01,  P-2  =  4.39,  P-3  =  5.79  were  obtained  using  initial  estimates  of  1. 

Linear  Versus  Nonlinear  Strain  Relations 

It  was  noted  that  for  k  =  2.7  N / /jrn  and  do  =  —10.6  MPa,  |^|^|  <<  | do |  which  permits  the 
nonlinear  strain  relation  (47)  to  be  approximated  by  the  explicit  relation  (48).  Here  we  quantify  the 
accuracy  of  (48)  as  a  function  of  the  spring  constant  and  prestress  level  for  the  considered  actuator 
configuration. 

We  first  define  the  relative  errors 

\\pl-pn\\2 

eP=  "  ||pjV||2  "  xl0°  .  e£="  ||gA,||2"  ><  100 
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where  PN,eN  and  PL,eL  respectively  denote  the  polarization  and  strain  values  given  by  (48)  and 
(47)  and  ||  •  || 2  is  the  L2  norm.  The  relative  errors  obtained  by  using  the  estimated  polarization  values 
summarized  in  Table  2  are  compiled  in  Table  3  for  various  combinations  of  k  and  <70 •  It  is  observed 
that  for  prestresses  less  than  -32  MPa  and  spring  constants  less  than  27  N//mr,  the  linear  relation 
provides  suitable  accuracy. 


Prestress 

k  =  0.27  N/fxm 

k  =  2.7  N/frni 

k  =  5.4  N/ fj,m 

k  =  27  N/^im 

-2  MPa 

(0.0049,  0.0059) 

(0.0489,  0.0584) 

(0.0969,  0.1157) 

(0.4462,  0.5438) 

-4  MPa 

(0.0047,  0.0059) 

(0.0464,  0.0589) 

(0.0919,  0.1171) 

(0.4223,  0.5489) 

-8  MPa 

(0.0042,  0.0061) 

(0.0415,  0.0603) 

(0.0820,  0.1193) 

(0.3774,  0.5602) 

-16  MPa 

(0.0034,  0.0063) 

(0.0331,  0.0625) 

(0.0654,  0.1238) 

(0.3013,  0.5856) 

-32  MPa 

(0.0032,  0.0072) 

(0.0321,  0.0719) 

(0.0638,  0.1434) 

(0.2960,  0.6794) 

Table  3:  Relative  error  (ep,e£)  between  the  linear  relation  (48)  and  nonlinear  relation  (47). 


6.2  Model  Validation  for  Rate-Dependent  Bulk  PZT 

Soft  ferroelectric  compounds  have  lower  coercive  fields,  more  hysteresis  loss,  and  larger  susceptibil¬ 
ities  and  piezoelectric  constants  than  hard  compounds.  Here  we  illustrate  the  performance  of  the 
model  for  characterizing  the  rate-dependence  of  the  soft  PZT  material  PIC151  which  has  the  compo¬ 
sition  Pb(Ni1/3  Sb2/3)03-PbTi03-PbZr03  ternary  phase  system  that  includes  about  2-3%  Pb(Ni1/3 
Sb2/3)03  near  the  morphotropic  phase  boundary  of  PZT.  We  compare  to  data  reported  in  [96],  that 
was  collected  at  rates  of  0.01  Hz,  0.1  Hz  and  1  Hz,  as  shown  in  Figure  3.  The  PZT  in  this  case  was 
a  5  x  5  x  15  mm3  cube  that  was  cut  from  a  bulk  ferroelectric  speciman  and  electroded  at  the  ends. 
Hence  the  construction  and  electrode  configuration  differ  from  the  multilayer  actuator  discussed  in 
Section  6.1. 

It  is  observed  that  coercive  fields  increase  as  the  frequency  increases  and  that  both  the  polarization 
and  strain  continue  to  increase  when  the  field  is  reversed  at  its  maximal  values  of  2  kV/rnm.  This 
is  due  to  the  fact  that  dipoles  do  not  switch  instantaneously  as  modeled  by  the  relations  —  and 

jr  j  j  Tgo 

in  the  likelihood  relations  (34).  Hence  some  dipoles  will  continue  to  switch  from  —180°  to  90° 
and  90°  to  180°  during  the  initial  unloading  process.  This  rate-dependence  is  due  to  the  time  scale 
difference  between  dipole  kinetics  and  electrical  or  mechanical  loading  rates. 

The  initial  and  optimized  parameters  estimated  through  a  fit  to  0.1  Hz  and  1.0  Hz  loading  rates 
are  reported  in  Table  4  and  the  corresponding  model  fits  are  shown  in  Figure  25  along  with  the 
model  predictions  at  0.1  Hz.  It  is  observed  that  the  coercive  field  increase  and  rate-dependence  of 

the  polarization  is  accurately  modeled  at  all  three  frequencies.  The  modeled  polarization  and  strain 

behavior  are  also  quite  accurate  including  the  nearly  parabolic  and  anhysteretic  1  Hz  P-e  behavior 
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Table  4:  Initial  and  optimized  PZT  parameter  values  estimated  through  a  fit  to  the  0.1  Hz  and  1.0  Hz 
data  from  [96].  Density  parameters:  Fc  =  0.66  x  106,o/  =  0.5  x  106,crc  =  0.35  and  aq  =  1.02,  ao  = 
1.04,  a_i  =  0.98,  a_2  =  1.13,  a_3  =  1.23,  ft  =  0.77,  ft  =  0.86,  ft-i  =  1.34,  ft_2  =  1.80,  ft-3  = 
1.90,  P-4  =  2.36. 
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Figure  25:  Fit  of  the  polarization  and  strain  model  to  data  from  [96]  for  loading  rates  of  (d)-(f)  0.1  Hz 
and  (g)-(i)  1  Hz.  Model  prediction  for  (a)-(c)  0.01  Hz:  ( - )  experimental  data,  ( - )  model  fit. 


which  we  note  was  modeled  using  linear  domain-level  constitutive  relations  in  combination  with 
grain-level  switching  mechanisms.  It  is  observed  that  whereas  the  model  fit  is  quite  accurate  at  1  Hz, 
the  model  does  not  incorporate  the  degree  of  delayed  dipole  switching  as  is  exhibited  by  the  data. 
This  indicates  that  aspects  of  the  90°  switching  model  can  still  be  improved. 

A  comparison  between  the  quasistatic  hysteresis  and  butterfly  loops  for  the  prestressed  multi¬ 
layer  actuator  and  bulk  sample,  respectively  plotted  in  Figures  23  and  25,  reveal  differences  in  the 
remanence  values,  coercive  fields,  and  post-switching  slopes.  This  is  substantiated  by  differences  in 
the  optimized  parameter  values  in  Tables  2  and  4.  This  is  due  in  part  to  differences  in  the  actuator 
construction  and  drive  configuration  and  it  illustrates  the  necessity  of  calibrating  models  for  a  specific 
actuator  configuration. 

6.3  Model  Validation  for  PLZT 

As  noted  in  Section  2.2,  8/65/35  PLZT  has  a  rhombohedral  crystal  structure  so  switches  of  70.5° 
and  109.5°  produce  strains  and  changes  in  the  polarization.  This  can  be  modeled  two  ways  using 
the  homogenized  energy  model.  The  more  rigorous  is  to  modify  the  energy  landscape  so  that  it 
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Figure  26:  Fit  of  the  D-field  and  strain  model  to  data  from  Lynch  [61]:  (a)-(c)  do  =  —3  MPa 
and  (j)-(Z')  cro  =  —15  MPa.  Model  predictions  for  prestresses  of  (d)-(f)  do  =  —6  MPa  and  (g)-(i) 
(Tq  =  —10  MPa;  ( - )  experimental  data,  ( - )  model  fit  and  predictions. 


has  minima  at  a  =  ±180,  70.5  and  109.5.  Alternatively,  one  can  apply  the  tetragonal  model  with 
a  =  ±180  and  90  in  a  phenomenological  manner  to  model  the  stress-dependence  of  the  strains  and 
polarization.  We  illustrate  the  accuracy  of  the  latter  approach  through  fits  and  predictions  using 
PLZT  data  collected  at  various  prestresses  as  reported  in  [61]. 

To  estimate  model  parameters,  we  first  fit  the  polarization  and  strain  model  (41)  to  data  collected 
at  prestresses  of  do  =  —3  MPa  and  do  =  —15  MPa;  note  that  Lynch  reports  H-field  data  rather 
than  polarization  so  modeled  results  were  obtained  using  the  relation  D  =  £qE  ±  P.  This  yielded 
the  parameters  summarized  in  Table  5  and  the  model  fits  shown  in  Figure  26(a)-(c)  and  (j )-(•£).  The 
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model  with  these  parameters  was  then  used  to  predict  the  responses  for  prestresses  of  oo  =  ~6  MPa 
and  do  =  —10  MPa  yielding  the  results  shown  in  Figure  26(d)-(i).  The  model  underpredicts  the 
amount  of  non- 180°  switching  in  the  strain  model  with  do  =  —  3  MPa,  which  may  be  due  in  part  to 
the  use  of  tetragonal  rather  than  rhombohedral  energy  landscapes,  but  accurately  fits  and  predicts 
the  switching  behavior  at  other  stress  levels. 

It  is  observed  that  the  model  reasonably  characterizes  the  nearly  quadratic  D-e  behavior  despite 
the  fact  that  it  employs  linear  constitutive  relations  rather  than  quadratic  electrostrictive  relations. 
As  detailed  in  Section  2.2,  the  quadratic  strain  behavior  for  this  material  is  due  to  dipole  rotation 
rather  than  classical  electrostrictive  effects  so  the  homogenized  energy  model,  which  employs  energy 
landscapes  to  quantify  dipole  switching,  provides  reasonably  accurate  material  characterization. 
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Table  5:  Initial  and  optimized  PLZT  parameter  values  for  the  90°  model  fitting  two  intermediate 
prestress  levels  do  =  —  3MPa  and  do  =  —15  MPa. 


7  Concluding  Remarks 

Hysteresis,  which  is  intrinsic  to  ferroelectric  materials,  involves  multiple  spatial  and  temporal  scales. 
At  the  domain  level,  strains  or  changes  in  polarization  are  due  to  stress- induced  material  deformations 
or  field-induced  ion  movement  and  the  behavior  is  often  reversible  and  linear.  Field  or  stress-induced 
dipole  switching  at  the  grain  level  produces  irreversible  hysteresis  in  both  the  field-polarization 
and  field-strain  relations.  For  single  crystal  materials  comprised  of  a  single  grain,  the  switching  is 
typically  rapid  thus  producing  sharp  hysteresis  and  butterfly  loops  in  quasistatic  operating  regimes. 
For  poly  crystalline  materials  with  distributed  interaction  and  coercive  fields,  hysteresis  and  butterfly 
loops  are  smoothed  due  to  nonuniform  grain  contributions.  The  behavior  of  hysteresis  loops  is  further 
modified  when  hysteretic  actuator  or  sensor  materials  are  employed  on  distributed  structures.  Hence 
the  nature  of  the  hysteretic  response  is  highly  influenced  by  the  spatial  scale  under  consideration. 

Furthermore,  ferroelectric  materials  exhibit  creep  and  rate-dependent  effects  even  at  low  frequen¬ 
cies.  This  is  due  to  the  fact  that  the  kinetics  associated  with  dipole  switching  typically  differs  from 
mechanical  or  electrical  loading  rates.  This  establishes  multiple  time  scales  that  must  be  incorporated 
in  dynamic  models. 

The  homogenized  energy  model  is  a  multiscale,  microscopically-motivated  or  nricronrechanical 
approach  that  incorporates  the  rate-dependence  and  multiple  timescales  exhibited  by  materials.  At 
the  domain  level,  the  minimization  of  Gibbs  energy  densities  yields  linear  constitutive  relations. 
At  the  grain  level,  dipole  fractions  serve  as  appropriate,  physically-motivated,  internal  variables 
to  quantify  the  timescales  and  hysteresis  associated  with  dipole  switching.  The  dynamics  of  dipole 
fractions  are  governed  by  evolution  equations  driven  by  likelihood  rates  constructed  using  Boltzmann 
theory  to  quantify  the  scaled  probability  of  transitioning  between  stable  equilibria  associated  with 
dipole  variants.  Macroscale  models  are  developed  by  assuming  that  properties  such  as  coercive  fields, 
critical  driving  forces,  and  interaction  fields  are  manifestations  of  underlying  densities  rather  than 
constants.  Finally,  it  is  shown  that  the  homogenized  energy  model  framework  facilitates  subsequent 
development  of  distributed  system  models.  The  complete  multiscale  development  is  illustrated  in 
Figure  4. 
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In  Section  2,  significant  discussion  was  devoted  to  mechanisms  that  produce  the  quadratic  P-e 
or  D-e  behavior  exhibited  by  ferroelectric  materials.  For  compounds  such  as  PMN,  this  is  due  to 
classical  electrostriction  and  quadratic  domain-level  constitutive  relations  are  required  to  model  the 
nearly  anhysteretic  field-strain  behavior.  In  other  materials,  such  as  PZT,  this  is  due  to  domain 
rotation  and  the  model  validation  results  illustrate  that  highly  accurate  material  characterization 
can  be  achieved  using  linear  domain  level  relations. 

The  validation  results  of  Section  6  demonstrate  the  capability  of  the  model  to  characterize  ma¬ 
jor  and  biased  minor  polarization  and  strain  loops,  the  effect  of  prestresses,  and  creep  and  rate- 
dependencies  for  PZT  and  PLZT.  Additional  validation  results  demonstrating  the  performance  of 
the  framework  for  characterizing  single  crystal  BaTiC>3  and  variable  loading  rates  in  PZT  are  pre¬ 
sented  in  the  companion  paper  [31].  It  is  also  shown  in  [31]  that  due  to  its  energy  basis,  the  model 
admits  highly  efficient  implementation  and  data-driven  algorithms  to  determine  initial  parameter 
estimates  based  on  measured  properties  of  the  data.  This  facilitates  model  calibration  and  imple¬ 
mentation  for  design  and  control  of  devices  and  complex  structures  arising  in  applications. 

We  note  that  the  polarization  model  (26)  bears  some  resemblance  to  a  Preisach  model.  In  the 
absence  of  thermal  activation,  (26)  satisfies  congruency  and  deletion  properties  and  hence  it  pro¬ 
vides  an  energy-based  derivation  for  certain  classical  Preisach  models.  As  detailed  in  [75],  extended 
Preisach  formulations  have  been  developed  to  incorporate  rate  effects,  thermal  activation  and  creep. 
However,  they  do  so  by  modifying  the  densities  rather  than  using  energy  and  kinetics  principles  to 
incorporate  them  in  the  kernels  as  is  the  case  for  the  homogenized  energy  model.  Moreover,  the 
combined  strain  and  polarization  relations  that  result  from  the  inclusion  of  elastic,  electric,  and  cou¬ 
pling  components  in  the  domain-level  Gibbs  energy  density  yield  macroscale  constitutive  relations 
that  incorporate  significantly  more  physics  than  Preisach-based  approaches. 

An  issue  that  is  of  significant  interest  for  control  design  of  systems  with  hysteretic  actuators 
and  sensors  concerns  the  construction  of  approximate  inverse  models.  For  the  homogenized  energy 
framework,  initial  approximate  inverse  models  have  been  implemented  at  rates  that  are  proven  no 
slower  than  1/6-1/7  the  rate  of  forward  simulations.  The  completion  of  these  algorithms  and  their 
incorporation  in  robust  control  designs  constitutes  present  and  future  research. 

A  Available  Codes 

To  facilitate  model  validation  and  dissemination  to  the  community,  we  have  made  codes  and  data 
available  at  the  website  http://www4.ncsu.edu/~jhcrews/smart/code/pzt/index.html. 
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